text:
##################################################################
"""
$(SIGNATURES)

Extract value from dual number. Use to debug physics callbacks.
Re-exported from ForwardDiff.jl
"""
const value=ForwardDiff.value

##################################################################
"""
$(SIGNATURES)

Extract derivatives from dual number.
"""
const partials=ForwardDiff.partials

##################################################################
"""
$(SIGNATURES)

Extract number of derivatives from dual number.
"""
const npartials=ForwardDiff.npartials

"""
Add value to matrix if it is nonzero
"""
@inline function _addnz(matrix,i,j,v::Tv,fac) where Tv
    if isnan(v)
        error("trying to assemble NaN")
    end
    if v!=zero(Tv)
        rawupdateindex!(matrix,+,v*fac,i,j)
    end
end

"""
$(TYPEDEF)

Exception thrown if Newton's method convergence fails.
"""
struct ConvergenceError <: Exception
end

"""
$(TYPEDEF)

Exception thrown if error occured during assembly (e.g. domain error)
"""
struct AssemblyError <: Exception
end

"""
$(TYPEDEF)

Exception thrown if error occured during factorization.
"""
struct FactorizationError <: Exception
end

"""
$(TYPEDEF)

Exception thrown if embedding fails
"""
struct EmbeddingError <: Exception
end

"""
Print error when catching exceptions
"""
function _print_error(err,st)
    println()
    println(err)
    nlines=5
    for i=1:min(nlines,length(st))
        line=@sprintf("%s",st[i])
        L=length(line)
        if L<80
            println(line)
        else
            print(line[1:35])
            print(" ... ")
            println(line[L-35:L])
        end
    end
    if length(st)>nlines
        println("...")
    end
    println()
end

"""
Solve time step problem. This is the core routine
for implicit Euler and stationary solve
"""
function _solve!(
    solution::AbstractMatrix{Tv}, # old time step solution resp. initial value
    oldsol::AbstractMatrix{Tv}, # old time step solution resp. initial value
    system::AbstractSystem{Tv}, # Finite volume system
    control::NewtonControl,
    tstep::Tv,
    log::Bool
) where Tv

    _complete!(system, create_newtonvectors=true)
    nlhistory=zeros(0)

    solution.=oldsol
    residual=system.residual
    update=system.update
    _initialize!(solution,system)
    SuiteSparse.UMFPACK.umf_ctrl[3+1]=control.umfpack_pivot_tolerance
    if typeof(control.factorization)!=typeof(system.factorization)
        system.factorization=control.factorization
    end
    oldnorm=1.0
    converged=false
    if control.verbose
        @printf("    Start Newton iteration\n")
    end
    nlu=0
    nround=0
    damp=control.damp_initial
    tolx=0.0
    rnorm=LinearAlgebra.norm(values(solution),1)

    for ii=1:control.max_iterations
        try
            eval_and_assemble(system,solution,oldsol,residual,tstep,edge_cutoff=control.edge_cutoff)
        catch err
            if (control.handle_exceptions)
                _print_error(err,stacktrace(catch_backtrace()))
                throw(AssemblyError())
            else
                rethrow(err)
            end
        end

        mtx=system.matrix

        nliniter=0
        # Sparse LU factorization
        if nlu==0 # (re)factorize, if possible reusing old factorization data
            try
                factorize!(system.factorization,mtx)
            catch err
                if (control.handle_exceptions)
                    _print_error(err,stacktrace(catch_backtrace()))
                    throw(FactorizationError())
                else
                    rethrow(err)
                end
            end
        end
        if issolver(system.factorization) && nlu==0
            # Direct solution via LU solve
            try
                ldiv!(values(update),system.factorization,values(residual))
            catch err
                if (control.handle_exceptions)
                    _print_error(err,stacktrace(catch_backtrace()))
                    throw(FactorizationError())
                else
                    rethrow(err)
                end
            end
        elseif !issolver(system.factorization) || nlu>0
            # Iterative solution
            update.=zero(Tv)
            (sol,history)= bicgstabl!(values(update),
                                      mtx,
                                      values(residual),
                                      1,
                                      Pl=system.factorization,
                                      reltol=control.tol_linear,
                                      max_mv_products=100,
                                      log=true)
            nliniter=history.iters
        else
            error("This should not happen")
        end

        nlu=min(nlu+1,control.max_lureuse)
        solval=values(solution)
        solval.-=damp*values(update)
        damp=min(damp*control.damp_growth,1.0)
        norm=LinearAlgebra.norm(values(update),Inf)
        if tolx==0.0
            tolx=norm*control.tol_relative
        end
        dnorm=1.0
        rnorm_new=LinearAlgebra.norm(values(solution),1)
        if rnorm>1.0e-50
            dnorm=abs((rnorm-rnorm_new)/rnorm)
        end

        if dnorm<control.tol_round
            nround=nround+1
        else
            nround=0
        end

        push!(nlhistory,norm)
        if control.verbose
            if   control.tol_linear<1.0
                itstring=@sprintf("it=% 3d(% 2d)",ii,nliniter)
            else
                itstring=@sprintf("it=% 3d",ii)
            end
            if control.max_round>0
                @printf("    %s du=%.3e cont=%.3e dnorm=%.3e %d\n",itstring,norm, norm/oldnorm,dnorm,nround)
            else
                @printf("    %s du=%.3e cont=%.3e\n",itstring,norm, norm/oldnorm)
            end
        end
        if ii>1 &&  norm/oldnorm > 1.0/control.tol_mono
            converged=false
            break
        end

        if norm<control.tol_absolute || norm <tolx
            converged=true
            break
        end
        oldnorm=norm
        rnorm=rnorm_new

        if nround>control.max_round
            converged=true
            break
        end
    end
    if !converged
        throw(ConvergenceError())
    end
    if control.verbose
        @printf("    Newton iteration successful\n")
    end
    return nlhistory
end

################################################################
"""
$(SIGNATURES)

Main assembly method.

Evaluate solution with result in right hand side F and
assemble Jacobi matrix into system.matrix.
"""
function eval_and_assemble(system::AbstractSystem{Tv, Ti},
                           U::AbstractMatrix{Tv}, # Actual solution iteration
                           UOld::AbstractMatrix{Tv}, # Old timestep solution
                           F::AbstractMatrix{Tv},# Right hand side
                           tstep::Tv; # time step size. Inf means stationary solution
                           edge_cutoff=0.0
                           ) where {Tv, Ti}

    _complete!(system) # needed here as well for test function system which does not use newton

    grid    = system.grid
    physics = system.physics
    node    = Node{Tv,Ti}(system)
    bnode   = BNode{Tv,Ti}(system)
    edge    = Edge{Tv,Ti}(system)
    bedge   = BEdge{Tv,Ti}(system)

    @create_physics_wrappers(physics, node, bnode, edge, bedge)

    nspecies::Int = num_species(system)
    matrix        = system.matrix

    cellnodefactors::Array{Tv,2}  = system.cellnodefactors
    celledgefactors::Array{Tv,2}  = system.celledgefactors
    bfacenodefactors::Array{Tv,2} = system.bfacenodefactors
    bfaceedgefactors::Array{Tv,2} = system.bfaceedgefactors

    # Reset matrix + rhs
    nzv  = nonzeros(matrix)
    nzv .= 0.0
    F   .= 0.0

    # structs holding diff results for storage, reaction,  flux ...
    result_r    = DiffResults.DiffResult(Vector{Tv}(undef, nspecies), Matrix{Tv}(undef, nspecies, nspecies))
    result_s    = DiffResults.DiffResult(Vector{Tv}(undef, nspecies), Matrix{Tv}(undef, nspecies, nspecies))
    result_flx  = DiffResults.DiffResult(Vector{Tv}(undef, nspecies), Matrix{Tv}(undef, nspecies, 2*nspecies))
    result_bflx = DiffResults.DiffResult(Vector{Tv}(undef, nspecies), Matrix{Tv}(undef, nspecies, 2*nspecies))
    # Array holding function results
    Y = Array{Tv,1}(undef,nspecies)

    # Arrays for gathering solution data
    UK    = Array{Tv,1}(undef,nspecies)
    UKOld = Array{Tv,1}(undef,nspecies)
    UKL   = Array{Tv,1}(undef,2*nspecies)

    # array holding source term
    src   = zeros(Tv,nspecies)
    bsrc  = zeros(Tv,nspecies)
    # arrays holding storage terms for old solution
    oldstor   = zeros(Tv, nspecies)
    res_react = zeros(Tv, nspecies)
    jac_react = zeros(Tv, nspecies, nspecies)
    oldbstor  = zeros(Tv, nspecies)

    # Due to

    # https://github.com/JuliaDiff/ForwardDiff.jl/issues/516

    # in   order   to  avoid   allocations,   we   directly  call   into
    # vector_mode_jacobian!. However, by default,  this assumes that the
    # length    of     the    argument     vector    is     less    than
    # DEFAULT_CHUNK_THRESHOLD.

    # This threshold can be increased by passing ForwardDiff.Chunk(UK,...).

    # See also https://juliadiff.org/ForwardDiff.jl/stable/user/advanced/#Configuring-Chunk-Size

    cfg_r     = ForwardDiff.JacobianConfig(reactionwrap, Y, UK, ForwardDiff.Chunk(UK,nspecies))
    cfg_s     = ForwardDiff.JacobianConfig(storagewrap, Y, UK, ForwardDiff.Chunk(UK,nspecies))
    cfg_br    = ForwardDiff.JacobianConfig(breactionwrap, Y, UK, ForwardDiff.Chunk(UK,nspecies))
    cfg_bs    = ForwardDiff.JacobianConfig(bstoragewrap, Y, UK, ForwardDiff.Chunk(UK,nspecies))

    cfg_bflx  = ForwardDiff.JacobianConfig(bfluxwrap, Y, UKL,ForwardDiff.Chunk(UKL,2*nspecies))
    cfg_flx   = ForwardDiff.JacobianConfig(fluxwrap, Y, UKL,ForwardDiff.Chunk(UKL,2*nspecies))

    # Inverse of timestep
    # According to Julia documentation, 1/Inf=0 which
    # comes handy to write compact code here for the
    # case of stationary problems.
    tstepinv = 1.0/tstep

    boundary_factors::Array{Tv,2} = system.boundary_factors
    boundary_values::Array{Tv,2}  = system.boundary_values
    bfaceregions::Vector{Ti} = grid[BFaceRegions]

    nbfaces = num_bfaces(grid)
    ncells  = num_cells(grid)
    geom=grid[CellGeometries][1]
    bgeom   = grid[BFaceGeometries][1]

    nn::Int = num_nodes(geom)
    ne::Int = num_edges(geom)

    ncalloc=@allocated  for icell=1:ncells
        for inode=1:nn
            _fill!(node,inode,icell)
            @views UK    .= U[:,node.index]
            @views UKOld .= UOld[:,node.index]
            # xx gather:
                # ii=0
                # for i region_spec.colptr[ireg]:region_spec.colptr[ireg+1]-1
                #    ispec=Fdof.rowval[idof]
                #    ii=ii+1
                #    node.speclist[ii]=ispec
                #    UK[ii]=U[ispec,K]
                #   UK[1:nspecies]=U[:,node.index]
                #   UKOld[1:nspecies]=UOld[:,node.index]
                # Evaluate source term

            if issource
                sourcewrap(src)
            end

            if isreaction
                # Evaluate & differentiate reaction term if present
                Y .= zero(Tv)
                ForwardDiff.vector_mode_jacobian!(result_r,reactionwrap,Y,UK,cfg_r)
                res_react = DiffResults.value(result_r)
                jac_react = DiffResults.jacobian(result_r)
            end

            # Evaluate & differentiate storage term
            Y.=zero(Tv)
            ForwardDiff.vector_mode_jacobian!(result_s,storagewrap,Y,UK,cfg_s)
            res_stor = DiffResults.value(result_s)
            jac_stor = DiffResults.jacobian(result_s)

            # Evaluate storage term for old timestep
            oldstor.=zero(Tv)
            storagewrap(oldstor,UKOld)
            K=node.index
            for idof=_firstnodedof(F,K):_lastnodedof(F,K)
                ispec=_spec(F,idof,K)
                _add(F,idof,cellnodefactors[inode,icell]*(res_react[ispec]-src[ispec] + (res_stor[ispec]-oldstor[ispec])*tstepinv))
                for jdof=_firstnodedof(F,K):_lastnodedof(F,K)
                    jspec=_spec(F,jdof,K)
                    _addnz(matrix,idof,jdof,jac_react[ispec,jspec]+ jac_stor[ispec,jspec]*tstepinv,cellnodefactors[inode,icell])
                end
            end
        end

        for iedge=1:ne
            if abs(celledgefactors[iedge,icell])<edge_cutoff
                continue
            end
            _fill!(edge,iedge,icell)

            #Set up argument for fluxwrap
            @views UKL[1:nspecies]            .= U[:, edge.node[1]]
            @views UKL[nspecies+1:2*nspecies] .= U[:, edge.node[2]]

            Y.=zero(Tv)
            ForwardDiff.vector_mode_jacobian!(result_flx,fluxwrap,Y,UKL,cfg_flx)
            res=DiffResults.value(result_flx)
            jac=DiffResults.jacobian(result_flx)

            K=edge.node[1]
            L=edge.node[2]

            fac=celledgefactors[iedge,icell]

            for idofK=_firstnodedof(F,K):_lastnodedof(F,K)
                ispec=_spec(F,idofK,K)
                idofL=dof(F,ispec,L)
                if idofL==0
                    continue
                end
                _add(F,idofK,fac*res[ispec])
                _add(F,idofL,-fac*res[ispec])

                for jdofK=_firstnodedof(F,K):_lastnodedof(F,K)
                    jspec=_spec(F,jdofK,K)
                    jdofL=dof(F,jspec,L)
                    if jdofL==0
                        continue
                    end

                    _addnz(matrix,idofK,jdofK,+jac[ispec,jspec         ],fac)
                    _addnz(matrix,idofL,jdofK,-jac[ispec,jspec         ],fac)
                    _addnz(matrix,idofK,jdofL,+jac[ispec,jspec+nspecies],fac)
                    _addnz(matrix,idofL,jdofL,-jac[ispec,jspec+nspecies],fac)
                 end
            end
        end
    end

    nbn::Int = num_nodes(bgeom)
    nbe::Int = num_edges(bgeom)

    # Assembly loop for boundary conditions
    nballoc= @allocated for ibface=1:nbfaces
        ibreg=bfaceregions[ibface]

        # Loop over nodes of boundary face
        for ibnode=1:nbn
            # Fill bnode data shuttle with data from grid
            _fill!(bnode,ibnode,ibface)

            # Copy unknown values from solution into dense array
            @views UK.=U[:,bnode.index]

            # Measure of boundary face part assembled to node
            bnode_factor::Tv=bfacenodefactors[ibnode,ibface]

            if isbsource
                bsourcewrap(bsrc)
            end

            # Global index of node
            K::Int=bnode.index

            # Assemble "standard" boundary conditions: Robin or
            # Dirichlet
            # valid only for interior species, currently not checked
            for ispec=1:nspecies
                idof=dof(F,ispec,K)
                # If species is present, assemble the boundary condition
                if idof>0
                    # Get user specified data
                    boundary_factor = boundary_factors[ispec,ibreg]
                    boundary_value  = boundary_values[ispec,ibreg]

                    if boundary_factor==Dirichlet
                        # Dirichlet is encoded in the boundary condition factor
                        # Penalty method: scale   (u-boundary_value) by penalty=boundary_factor

                        # Add penalty*boundry_value to right hand side
                        F[ispec,K]+=boundary_factor*(U[ispec,K]-boundary_value)

                        # Add penalty to matrix main diagonal (without bnode factor, so penalty
                        # is independent of h)
                        _addnz(matrix,idof,idof,boundary_factor,1)
                    else
                        # Robin boundary condition
                        F[ispec,K]+=bnode_factor*(boundary_factor*U[ispec,K]-boundary_value)
                        _addnz(matrix,idof,idof,boundary_factor, bnode_factor)
                    end
                end
            end

            # Boundary reaction term has been given.
            # Valid for both boundary and interior species
            if isbreaction
                # Evaluate function + derivative
                Y.=zero(Tv)
                ForwardDiff.vector_mode_jacobian!(result_r,breactionwrap,Y,UK,cfg_br)
                res_breact=DiffResults.value(result_r)
                jac_breact=DiffResults.jacobian(result_r)

                # Assemble RHS + matrix
                for idof=_firstnodedof(F,K):_lastnodedof(F,K)
                    ispec=_spec(F,idof,K)
                    _add(F,idof,bnode_factor*(res_breact[ispec]-bsrc[ispec]))
                    for jdof=_firstnodedof(F,K):_lastnodedof(F,K)
                        jspec=_spec(F,jdof,K)
                        _addnz(matrix,idof,jdof,jac_breact[ispec,jspec],bnode_factor)
                    end
                end

            end

            # Boundary reaction term has been given.
            # Valid only for boundary species, but system is currently not checked.
            if isbstorage

                # Fetch data for old timestep
                @views UKOld.=UOld[:,bnode.index]

                # Evaluate storage term for old timestep
                bstoragewrap(oldbstor,UKOld)

                # Evaluate & differentiate storage term for new time step value
                Y.=zero(Tv)
                ForwardDiff.vector_mode_jacobian!(result_s,bstoragewrap,Y,UK,cfg_bs)
                res_bstor=DiffResults.value(result_s)
                jac_bstor=DiffResults.jacobian(result_s)

                for idof=_firstnodedof(F,K):_lastnodedof(F,K)
                    ispec=_spec(F,idof,K)
                    # Assemble finite difference in time for right hand side
                    _add(F,idof,bnode_factor*(res_bstor[ispec]-oldbstor[ispec])*tstepinv)
                    # Assemble matrix.
                    for jdof=_firstnodedof(F,K):_lastnodedof(F,K)
                        jspec=_spec(F,jdof,K)
                        _addnz(matrix,idof,jdof,jac_bstor[ispec,jspec],bnode_factor*tstepinv)
                    end
                end
            end # if isbstorage

        end # ibnode=1:nbn
        if isbflux
            for ibedge=1:nbe
                if abs(bfaceedgefactors[ibedge,ibface]) < edge_cutoff
                    continue
                end

                _fill!(bedge,ibedge,ibface)
                @views UKL[1:nspecies]            .= U[:, bedge.node[1]]
                @views UKL[nspecies+1:2*nspecies] .= U[:, bedge.node[2]]

                Y.=zero(Tv)
                ForwardDiff.vector_mode_jacobian!(result_bflx, bfluxwrap, Y, UKL, cfg_bflx)
                res = DiffResults.value(result_bflx)
                jac = DiffResults.jacobian(result_bflx)

                K   = bedge.node[1]
                L   = bedge.node[2]

                fac = bfaceedgefactors[ibedge, ibface]
                #@show fac
                for idofK = _firstnodedof(F, K):_lastnodedof(F, K)
                    ispec =_spec(F, idofK, K)
                    idofL = dof(F, ispec, L)
                    if idofL == 0
                        continue
                    end
                    _add(F, idofK,  fac*res[ispec])
                    _add(F, idofL, -fac*res[ispec])

                    for jdofK = _firstnodedof(F,K):_lastnodedof(F,K)
                        jspec = _spec(F,jdofK,K)
                        jdofL = dof(F,jspec,L)
                        if jdofL == 0
                            continue
                        end

                        _addnz(matrix, idofK, jdofK, +jac[ispec, jspec         ], fac)
                        _addnz(matrix, idofL, jdofK, -jac[ispec, jspec         ], fac)
                        _addnz(matrix, idofK, jdofL, +jac[ispec, jspec+nspecies], fac)
                        _addnz(matrix, idofL, jdofL, -jac[ispec, jspec+nspecies], fac)
                     end
                end

            end
        end
    end

    # if  no new matrix entries have been created, we should see no allocations
    # in the previous two loops
    if isnothing(matrix.lnkmatrix) && !_check_allocs(system,ncalloc+nballoc)
        error("""Allocations in assembly loop: cells: $(ncalloc), bfaces: $(nballoc)
                            See the documentation of `check_allocs!` for more information""")
    end
   _eval_and_assemble_generic_operator(system,U,F)
   _eval_and_assemble_inactive_species(system,U,UOld,F)
end

"""
Evaluate and assemble jacobian for generic operator part.
"""
function _eval_and_assemble_generic_operator(system::AbstractSystem,U,F)

    if !has_generic_operator(system)
        return
    end
    generic_operator(f,u)=system.physics.generic_operator(f,u,system)
    vecF=values(F)
    vecU=values(U)
    y=similar(vecF)
    generic_operator(y,vecU)
    vecF.+=y
    forwarddiff_color_jacobian!(system.generic_matrix, generic_operator, vecU, colorvec = system.generic_matrix_colors)
    rowval=system.generic_matrix.rowval
    colptr=system.generic_matrix.colptr
    nzval=system.generic_matrix.nzval
    for i=1:length(colptr)-1
        for j=colptr[i]:colptr[i+1]-1
            updateindex!(system.matrix,+,nzval[j],i,rowval[j])
        end
    end
end

################################################################
"""
````
solve!(solution, inival, system;
    control=NewtonControl(),
    tstep=Inf,
    log=false)
````

Perform solution of stationary problem(if `tstep==Inf`) or of  one step
of the implicit Euler method using Newton's method with `inival` as initial
value. The method writes into `solution`.

It returns `solution` or, if `log==true`, a tuple of solution and a vector
containing the residual history of Newton's method.
"""
function solve!(
    solution::AbstractMatrix{Tv}, # Solution
    inival::AbstractMatrix{Tv},   # Initial value
    system::AbstractSystem{Tv};     # Finite volume system
    control=NewtonControl(),      # Newton solver control information
    tstep::Tv=Inf,                # Time step size. Inf means  stationary solution
    log::Bool=false
) where Tv
    if control.verbose
        @time begin
            history=_solve!(solution,inival,system,control,tstep,log)
        end
    else
        history=_solve!(solution,inival,system,control,tstep,log)
    end
    return log ? (solution,history) : solution
end

################################################################
"""
````
solve(inival, system;
      control=NewtonControl(),
      tstep=Inf,
      log=false)
````

Perform solution of stationary problem(if `tstep==Inf`) or one step
of the implicit Euler method using Newton's method with `inival` as initial
value.
It returns a solution array or, if `log==true`, a tuple of solution and a vector
containing the residual history of Newton's method.
"""
function solve(
    inival::AbstractMatrix{Tv},   # Initial value
    system::AbstractSystem{Tv};     # Finite volume system
    control=NewtonControl(),      # Newton solver control information
    tstep::Tv=Inf,                # Time step size. Inf means  stationary solution
    log::Bool=false
) where Tv
    solve!(unknowns(system),inival,system,control=control, tstep=tstep,log=log)
end

################################################################
"""
````
function embed!(solution, inival, system;
                control=NewtonControl(),
                pre=function(sol,p) end,
                post=function(sol,p) end
)
````
Solve stationary problem via parameter embedding, calling
`solve!` for increasing values of the parameter p from interval
(0,1). The user is responsible for the interpretation of
the parameter.

The optional `pre()` callback can be used
to communicate its value to the system.

The optional `post()` callback method can be used to perform
various postprocessing steps.

If `control.handle_error` is true, if `solve!`  throws an error,
stepsize `control.Δp` is lowered, and `solve!` is called again with a smaller  parameter
value. If `control.Δp<control.Δp_min`, `embed!` is aborted
with error.

"""
function embed!(
    solution::AbstractMatrix{Tv}, # Solution
    inival::AbstractMatrix{Tv},   # Initial value
    system::AbstractSystem{Tv};     # Finite volume system
    control=NewtonControl(),      # Newton solver control information
    pre=function(sol,p) end,       # Function for preparing step
    post=function(sol,p) end      # Function for postprocessing successful step
) where Tv
    inival=copy(inival)
    p=0.0
    Δp=control.Δp
    if control.verbose
        @printf("  Embedding: start\n")
    end
    istep=0
    while p<1.0
        solved=false
        p0=p
        while !solved
            solved=true
            try
                p=min(p0+Δp,1.0)
                pre(solution,p)
                solve!(solution,inival, system ,control=control)
            catch err
                if (control.handle_exceptions)
                    _print_error(err,stacktrace(catch_backtrace()))
                else
                    rethrow(err)
                end
                # here we reduce the embedding step and retry the solution
                Δp=Δp*0.5
                if Δp<control.Δp_min
                    throw(EmbeddingError())
                end
                solved=false
                if control.verbose
                    @printf("  Embedding: retry: Δp=%.3e\n",Δp)
                end
            end
        end
        istep=istep+1
        if control.verbose
            @printf("  Embedding: istep=%d p=%.3e\n",istep,p)
        end
        Δp=min(control.Δp_max,Δp*1.2)
        post(solution,p)
        inival.=solution
    end
    if control.verbose
        @printf("  Embedding: success\n")
    end
    return solution
end

################################################################
"""
````
function evolve!(solution, inival, system, times;
                 control=NewtonControl(),
                 pre=function(sol,t) end,
                 post=function(sol,oldsol, t, Δt) end,
                 sample=function(sol,t) end,
                 delta=(u,v,t, Δt)->norm(sys,u-v,Inf)
)
````

Use implicit Euler method  + damped   Newton's method  to
solve time dependent problem. Time step control is performed
according to the data in `control`.  All times in `times`
are reached exactly.

Callbacks:
- `pre` is invoked before each time step
- `post`  is invoked after each time step
- `sample` is called for all times in `times[2:end]`.

`delta` is  used to control the time step.
"""
function evolve!(
    solution::AbstractMatrix{Tv}, # Solution
    inival::AbstractMatrix{Tv},   # Initial value
    system::AbstractSystem{Tv},    # Finite volume system
    times::AbstractVector;
    control=NewtonControl(),      # Newton solver control information
    pre=function(sol,t) end,       # Function for preparing step
    post=function(sol,oldsol, t, Δt) end,      # Function for postprocessing successful step
    sample=function(sol,t) end,      # Function to be called for each t\in times[2:end]
    delta=(u,v,t, Δt)->norm(sys,u-v,Inf) # Time step error estimator
) where Tv
    inival=copy(inival)
    _initialize_dirichlet!(inival,system)
    Δt=control.Δt
    if control.verbose
        @printf("  Evolution: start\n")
    end
    for i=1:length(times)-1

        Δt=max(Δt,control.Δt_min)
        tstart=times[i]
        tend=times[i+1]
        t=tstart
        istep=0

        while t<tend
            solved=false
            t0=t
            Δu=0.0
            while !solved
                solved=true
                try
                    t=t0+Δt
                    pre(solution,t)
                    solve!(solution,inival, system ,control=control,tstep=Δt)
                catch err
                    if (control.handle_exceptions)
                        _print_error(err,stacktrace(catch_backtrace()))
                    else
                        rethrow(err)
                    end
                    solved=false
                end
                if solved
                    Δu=delta(solution, inival,t, Δt)
                    if Δu>2.0*control.Δu_opt
                        solved=false
                    end
                end
                if !solved
                    # reduce time step and retry  solution
                    Δt=Δt*0.5
                    if Δt<control.Δt_min
                        @printf(" Δt_min=%.2g reached while Δu=%.2g >>  Δu_opt=%.2g\n",control.Δt_min, Δu,control.Δu_opt)
                        if !(control.force_first_step && istep==0)
                            throw(EmbeddingError())
                        else
                            solved=true
                        end
                    end
                    if control.verbose
                        @printf("  Evolution: retry: Δt=%.3e\n",Δt)
                    end
                end
            end
            istep=istep+1
            if control.verbose
                @printf("  Evolution: istep=%d t=%.3e Δu=%.3e\n",istep,t,Δu)
            end
            post(solution,inival,t, Δt)
            inival.=solution
            if t<tend
                Δt=min(control.Δt_max,
                       Δt*control.Δt_grow,
                       Δt*control.Δu_opt/(Δu+1.0e-14),
                       tend-t)
            end
        end
        sample(solution,times[i+1])
    end
    if control.verbose
        @printf("  Evolution: success\n")
    end
    return solution
end

"""
````
function solve(inival, system, times;
               control=NewtonControl(),
               pre=function(sol,t) end,
               post=function(sol,oldsol, t, Δt) end,
               sample=function(sol,t) end,
               delta=(u,v,t, Δt)->norm(sys,u-v,Inf),
               store_all=true,
               in_memory=true
)
````
Use implicit Euler method  + damped   Newton's method  to
solve time dependent problem. Time step control is performed
according to the data in `control`.  All times in the vector
`times` are reached exactly.

Callbacks:
- `pre` is invoked before each time step
- `post`  is invoked after each time step
- `sample` is called for all times in `times[2:end]`.

`delta` is  used to control the time step.

If `store_all==true`, all timestep solutions are stored. Otherwise,
only solutions for the moments defined in `times` are stored.

Returns a transient solution object `sol` containing the stored solution,
see [`TransientSolution`](@ref)

"""
function solve(inival::AbstractMatrix,
               sys::AbstractSystem,
               times::AbstractVector;
               control=NewtonControl(),
               pre=function(sol,t) end,       # Function for preparing step
               post=function(sol,oldsol, t, Δt) end,      # Function for postprocessing successful step
               sample=function(sol,t) end,      # Function to be called for each t\in times[2:end]
               delta=(u,v,t, Δt)->norm(sys,u-v,Inf), # Time step error estimator
               store_all=true, # if true, store all solutions, otherwise, store only at sampling times
               in_memory=true
               )
    tsol=TransientSolution(Float64(times[1]),inival, in_memory=in_memory)
    solution=copy(inival)
    if store_all
        evolve!(solution,inival,sys,times,
                control=control,
                post=(u,uold,t,Δt)->(post(u,uold,t,Δt);append!(tsol,t,u)),
                pre=pre,
                sample=sample,
                delta=delta
                )
    else
        evolve!(solution,inival,sys,times,
                control=control,
                post=post,
                sample=(u,t)->(sample(u,t);append!(tsol,t,u)),
                pre=pre,
                delta=delta
                )
    end
    tsol
end


text:
#! /usr/bin/env julia

f = "./test.txt"

function cb(a,b,c)
    @show a, b, c
    if b.changed
        println("$a changed")
    end
end

# I am not sure how to get this to work without calling watch_file again
while true
    a = watch_file(cb, f)
    wait(a)
end

ext: jl

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground


text:
function solve2(items, boxes)
    ret = 0
    for (w, v) in items
        bb = 10e6+1
        vv = -1
        ii = -1
        for (i, b) in enumerate(boxes)
            if b ≥ w && b ≤ bb
                bb = b
                vv = v
                ii = i
            end
        end
        if vv > 0
            ret += vv
            deleteat!(boxes, ii)
        end
    end
    ret
end

function solve()
    n, m, q = [parse(Int, x) for x in split(readline())]
    items = Array{Tuple{Int, Int}}([])
    for _ in 1:n
        item = Tuple(parse(Int, x) for x in split(readline()))
        push!(items, item)
    end
    sort!(items, by=x->(x[2], x[1]), rev=true)
    boxes = [parse(Int, x) for x in split(readline())]

    for _ in 1:q
        l, r = [parse(Int, x) for x in split(readline())]
        bs = append!(boxes[1:l-1], boxes[r+1:end])
        ret = solve2(items, bs)
        println(ret)
    end
end

solve()

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
function Base.sort!(idx::AbstractArray, vals::AbstractArray...)
    if !issorted(idx)
        ordering = sortperm(idx)
        permute!(idx, ordering)
        for v in vals
            permute!(v, ordering)
        end
    end
end

backed_matrix(obj::Union{HDF5.File, HDF5.Group}) = SparseDataset(obj)
backed_matrix(obj::HDF5.Dataset) = TransposedDataset(obj)

function hdf5_object_name(obj::Union{HDF5.File, HDF5.Group, HDF5.Dataset})
    name = HDF5.name(obj)
    return name[(last(findlast("/", name)) + 1):end]
end

function find_unique_colnames(mdata::MuData, property::Symbol, ncols::Int)
    nchars = 16
    allunique = false
    local colnames::Vector{String}
    while !allunique
        colnames = [randstring(nchars) for _ in 1:ncols]
        allunique = length(Set(colnames)) == ncols
        nchars *= 2
    end

    it = Iterators.flatten(((mdata,), values(mdata.mod)))
    for i in 1:ncols
        finished = false
        while !finished
            for ad in it
                try
                    names(getproperty(ad, property), colnames[i])
                    colnames[i] = "_" * colnames[i]
                    break
                catch e
                    if !(e isa ArgumentError)
                        rethrow(e)
                    end
                end
            end
            finished = true
        end
    end
    return colnames
end

function index_duplicates(idx::AbstractArray{T}) where T
    counter = Dict{T, UInt8}()
    dup_idx = Vector{UInt8}(undef, length(idx))
    for (i, val) in enumerate(idx)
        if val in keys(counter)
            count = counter[val] + 1
        else
            count = 0
        end
        dup_idx[i] = counter[val] = count
    end
    return dup_idx
end

function minimum_unsigned_type_for_n(n::Number)
    local mintype::Type
    for type in [UInt8, UInt16, UInt32, UInt64, UInt128]
        mval = typemax(type)
        if mval >= n
            mintype = type
            break
        end
    end
    return mintype
end

@inline function convertidx(
    idx::Union{AbstractUnitRange, Colon, AbstractVector{<:Integer}, AbstractVector{Bool}},
    ref::AbstractIndex{<:AbstractString},
)
    return idx
end
@inline function convertidx(idx::Number, ref::AbstractIndex{<:AbstractString})
    return idx:idx
end
@inline convertidx(
    idx::Union{AbstractString, AbstractVector{<:AbstractString}},
    ref::AbstractIndex{<:AbstractString},
) = ref[idx, true]

Base.axes(A::Union{AbstractMuData, AbstractAnnData}) = map(n -> Base.OneTo(n), size(A))
isbacked(ad::Union{AbstractMuData, AbstractAnnData}) = !isnothing(file(ad))

@inline function Base.checkbounds(::Type{Bool}, A::Union{AbstractMuData, AbstractAnnData}, I...)
    Base.checkbounds_indices(
        Bool,
        axes(A),
        Tuple(i isa AbstractString || i isa AbstractVector{<:AbstractString} ? (:) : i for i in I),
    )
end

@inline function Base.checkbounds(A::Union{AbstractMuData, AbstractAnnData}, I...)
    checkbounds(Bool, A, I...) || throw(BoundsError(A, I))
    nothing
end

function Base.summary(A::Union{AbstractMuData, AbstractAnnData})
    s = size(A)
    return "$(typeof(A)) with $(s[1]) observations and $(s[2]) variables"
end

Base.summary(io::IO, A::Union{AbstractMuData, AbstractAnnData}) = print(io, summary(A))

Base.firstindex(A::Union{AbstractMuData, AbstractAnnData}, d::Integer) = 1
Base.lastindex(A::Union{AbstractMuData, AbstractAnnData}, d::Integer) = size(A, d)

Base.copy(d::Union{MuDataView, AnnDataView}) = parent(d)[parentindices(d)...]


text:
import Base:in, +, -, *, ==, convert, isempty, isequal, length, push!, sizeof, append!
import SeisIO: GeoLoc, PZResp
export RawData

const rawfields = [:name, :loc, :fs, :gain, :freqmin, :freqmax, :cc_step,
                   :cc_len, :time_norm, :resp]

"""
   RawData

A structure for raw ambient noise data.

## Fields: RawData
| **Field** | **Description** |
|:--------------------|:--------------------|
| :name       | Freeform channel names |
| :id         | Channel ids. use NET.STA.LOC.CHAN format when possible. |
| :loc        | Location (position) object |
| :fs         | Sampling frequency in Hz. |
| :gain       | Scalar gain; divide data by the gain to convert to units  |
| :freqmin    | Minimum frequency from instrument response.  |
| :freqmax    | Maximum frequency from instrument response. |
| :cc_len     | Raw data window length in seconds. |
| :cc_step    | Spacing between windows in seconds. |
| :time_norm  | Apply one-bit whitening with "one_bit". |
| :resp       | Instrument response object, format [zeros poles] |
| :misc       | Dictionary for non-critical information. |
| :notes      | Timestamped notes; includes automatically-logged acquisition and |
|             | processing information. |
| :t          | Starttime of each window. |
| :x          | Raw data stored in columns. |
"""
mutable struct RawData
  name::String                                # name [Net.Sta.Loc.Chan]
  id::String                                  # id [Y-mm-dd] this is date of raw
  loc::GeoLoc                                 # loc
  fs::Float64                                 # sampling rate [Hz]
  gain::Float64                               # gain
  freqmin::Float64                            # minumum frequency [Hz]
  freqmax::Float64                            # maximum frequency [Hz]
  cc_len::Int                                 # window_length [s]
  cc_step::Int                                # step between windows [s]
  whitened::Bool                              # whitening applied
  time_norm::Union{Bool,String}               # time normaliation
  resp::PZResp                                # response poles/zeros
  misc::Dict{String,Any}                      # misc
  notes::Array{String,1}                      # notes
  t::Array{Float64,1}                         # time
  x::Array{<:Union{Float32,Float64},2}  # raw data

  function RawData(
      name     ::String,
      id       ::String,
      loc      ::GeoLoc,
      fs       ::Float64,
      gain     ::Float64,
      freqmin  ::Float64,
      freqmax  ::Float64,
      cc_len   ::Int,
      cc_step  ::Int,
      whitened ::Bool,
      time_norm::Union{Bool,String},
      resp     ::PZResp,
      misc     ::Dict{String,Any},
      notes    ::Array{String,1},
      t        ::Array{Float64,1},
      x        ::Array{<:Union{Float32,Float64},2}
      )

      return new(name, id, loc, fs, gain, freqmin, freqmax, cc_len, cc_step,
                 whitened, time_norm, resp, misc, notes, t, x)
    end

    function RawData(S::SeisData,cc_len::Int,cc_step::Int)
      merge!(S)
      # subset by time
      starttime, endtime = u2d.(nearest_start_end(S[1],cc_len, cc_step))
        # check if waveform length is < cc_len
      if Int(floor((endtime - starttime).value / 1000)) < cc_len
          return nothing
      end
      # sync start and end times
      x, starts = slide(sync(S[1],s=starttime,t=endtime) , cc_len, cc_step)
      return new(S[1].id,Dates.format(u2d(S[1].t[1,2]*1e-6),"Y-mm-dd"),S[1].loc,
                 S[1].fs,S[1].gain,1/cc_len,S[1].fs/2,cc_len,cc_step,false,
                 false,S[1].resp,S[1].misc,S[1].notes,starts,x)
    end

    function RawData(C::SeisChannel,cc_len::Int,cc_step::Int)
      # subset by time
      starttime, endtime = u2d.(nearest_start_end(C,cc_len, cc_step))
        # check if waveform length is < cc_len
      if Int(floor((endtime - starttime).value / 1000)) < cc_len
          return nothing
      end
      # sync start and end times
      x, starts = slide(sync(C,s=starttime,t=endtime), cc_len, cc_step)
      return new(C.id,Dates.format(u2d(C.t[1,2]*1e-6),"Y-mm-dd"),C.loc,C.fs,
                 C.gain,1/cc_len,C.fs/2,cc_len,cc_step,false,false, C.resp,
                 C.misc,C.notes,starts,x)
    end
end

RawData(;
          name     ::String                    = "",
          id       ::String                    = "",
          loc      ::GeoLoc                    = GeoLoc(),
          fs       ::Float64                   = zero(Float64),
          gain     ::Float64                   = one(Float64),
          freqmin  ::Float64                   = zero(Float64),
          freqmax  ::Float64                   = zero(Float64),
          cc_len   ::Int                       = zero(Int),
          cc_step  ::Int                       = zero(Int),
          whitened ::Bool                      = false,
          time_norm::Union{Bool,String}        = false,
          resp     ::PZResp                    = PZResp(),
          misc     ::Dict{String,Any}          = Dict{String,Any}(),
          notes    ::Array{String,1}           = Array{String,1}(undef, 0),
          t        ::Array{Float64,1}          = Array{Float64,1}(undef, 0),
          x      ::Array{<:Union{Float32,Float64},2} =
                     Array{Float32,2}(undef, 0, 2)
          ) = RawData(name, id, loc, fs, gain, freqmin, freqmax, cc_len, cc_step,
                      whitened, time_norm, resp, misc, notes, t, x)

in(s::String, R::RawData) = R.id==s

isempty(R::RawData) = min(R.name=="",R.id == "",isempty(R.loc),
                          R.fs == zero(typeof(F.fs)),
                          R.gain == one(typeof(F.gain)),
                          R.freqmin == zero(typeof(F.freqmin)),
                          R.freqmax == zero(typeof(F.freqmax)),
                          R.cc_len == zero(typeof(F.cc_len)),
                          R.cc_step == zero(typeof(F.cc_step)),
                          F.time_norm == false,
                          isempty(R.resp), isempty(R.misc),isempty(R.notes),
                          isempty(R.t),isempty(R.raw))

isequal(R::RawData, U::RawData) = minimum([hash(getfield(R,i))==hash(getfield(U,i)) for i in rawfields]::Array{Bool,1})
==(F::RawData, U::RawData) = isequal(R,U)::Bool

append!(R::RawData, U::RawData)  = (if R == U;
  setfield!(R, :t, vcat(getfield(R,:t), getfield(U,:t)));
  setfield!(R, :raw, hcat(getfield(R,:raw), getfield(U,:raw)))
  elseif isempty(R);
  [setfield!(R, i, getfield(U,i)) for i in rawfields];
  [setfield!(R, i, getfield(U,i)) for i in [:t,:raw]];
  end;
  return R)
+(R::RawData, U::RawData) = (T = deepcopy(R); return append!(T, U))

function sizeof(R::RawData)
s = sum([sizeof(getfield(R,f)) for f in rawfields])
if !isempty(R.notes)
  s += sum([sizeof(i) for i in R.notes])
end
if !isempty(R.misc)
  s += sum([sizeof(i) for i in values(R.misc)])
end
return s
end


text:
using Clang.JLLEnvs
using Clang.Generators
using JpegTurbo_jll

cd(@__DIR__)

include_dir = normpath(JpegTurbo_jll.artifact_dir, "include")

function generate_prologue(options)
    # Macro "#define LIBJPEG_TURBO_VERSION 2.1.0" is not well-defined values
    # for Clang to parse. Thus we manually parse and generate it as prologue.
    function get_jpegturbo_version()
        lines = readlines(joinpath(include_dir, "jconfig.h"))
        prefix = "#define LIBJPEG_TURBO_VERSION "
        version_line = filter(lines) do line
            startswith(line, prefix)
        end[1]
        version_str = strip(split(version_line, prefix)[2])
        return VersionNumber(version_str)
    end

    outfile = options["general"]["prologue_file_path"]
    open(outfile, "w") do io
        println(io, "const LIBJPEG_TURBO_VERSION = v\"$(get_jpegturbo_version())\"")
        println(io)

        # https://github.com/libjpeg-turbo/libjpeg-turbo/blob/14ce28a92d45e4e22b643bd845ba6c543ebcd388/win/jconfig.h.in#L14-L18
        println(io, "const boolean = @static Sys.iswindows() ? Cuchar : Cint")
        println(io)
    end
end

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$include_dir")
push!(args, "--include=stdio.h")

generate_prologue(options)

# https://github.com/JuliaInterop/Clang.jl/discussions/373
headers = [
    joinpath(include_dir, "jpeglib.h"),
    joinpath(include_dir, "turbojpeg.h")
]

ctx = create_context(headers, args, options)

build!(ctx)


text:
module FNVHash

export fnv1, fnv1a, fnv_offset_basis, fnv_prime

const FNVInt = Union{UInt32, UInt64, UInt128}

@inline fnv_offset_basis(::Type{UInt32}) = 0x811c9dc5
@inline fnv_offset_basis(::Type{UInt64}) = 0xcbf29ce484222325
@inline fnv_offset_basis(::Type{UInt128}) = 0x6c62272e07bb014262b821756295c58d

@inline fnv_prime(::Type{UInt32}) = 0x01000193
@inline fnv_prime(::Type{UInt64}) = 0x00000100000001B3
@inline fnv_prime(::Type{UInt128}) = 0x0000000001000000000000000000013B

"""
    fnv1(T, p::Ptr{UInt8}, n::Integer)
    fnv1(T, x::Union{String, DenseVector{UInt8}})

Calculate the Fowler-Noll-Vo version 1 hash of the input for the bitwidth defined by T. T must unsigned.
"""
function fnv1 end

"""
    fnv1a(T, p::Ptr{UInt8}, n::Integer)
    fnv1a(T, x::Union{String, DenseVector{UInt8}})

Calculate the Fowler-Noll-Vo version 1a hash of the input for the bitwidth defined by T. T must unsigned.
"""
function fnv1a end

function fnv1(::Type{T}, p::Ptr{UInt8}, n::Integer) where T <: FNVInt
    h = fnv_offset_basis(T)
    for i in 1:n
        x = unsafe_load(p, i)
        h = (h * fnv_prime(T)) ⊻ x
    end
    h
end

function fnv1a(::Type{T}, p::Ptr{UInt8}, n::Integer) where T <: FNVInt
    h = fnv_offset_basis(T)
    for i in 1:n
        x = unsafe_load(p, i)
        h = (h ⊻ x) * fnv_prime(T)
    end
    h
end

fnv1(::Type{T}, s::Union{String, DenseVector{UInt8}}) where T <: FNVInt = GC.@preserve s fnv1(T, pointer(s), sizeof(s))
fnv1a(::Type{T}, s::Union{String, DenseVector{UInt8}}) where T <: FNVInt = GC.@preserve s fnv1a(T, pointer(s), sizeof(s))

end


text:
using JuliaProject

@test 2 == double(1)
input = "x"
@test input * input == double(input)


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

# This code was ported from a MATLAB implementation available at
# http://www.piso.at/julius/index.php/projects/programmierung/13-2d-wave-equation-in-octave
# with permission of the author.
# Original MATLAB implementation is copyright (c) 2014 Julius Piso.

using ParallelAccelerator
using DocOpt
if "--demo" in ARGS
    using Winston
end

@acc function runWaveStencil(p, c, f, r, t, n, s2, s4, s)
    runStencil(p, c, f, 1, :oob_skip) do p, c, f
        f[0, 0] = 2 * c[0, 0] - p[0, 0] + r * r * (c[-1, 0] + c[1, 0] + c[0, -1] + c[0, 1] - 4*c[0, 0])
    end

    # Dynamic source
    f[s2+s4-1:s2+s4+1, 1:2] = 1.5 * sin.(t*n)
    f[s2-s4-1:s2-s4+1, 1:2] = 1.5 * sin.(t*n)
    f[2:s-1, 1:2] = 1.0

    # Transparent boundary handling
    f[2:s-1, 1:1] = (2.0 * c[2:s-1, 1:1] + (r-1.0) * p[2:s-1, 1:1] + 2.0*r*r*(c[2:s-1, 2:2] + c[3:s, 1:1] + c[1:s-2, 1:1] - 3.0 * c[2:s-1, 1:1])) / (1.0+r) # Y:1
    f[2:s-1, s:s] = (2.0 * c[2:s-1, s:s] + (r-1.0) * p[2:s-1, s:s] + 2.0*r*r*(c[2:s-1, s-1:s-1] + c[3:s, s:s] + c[1:s-2, s:s] - 3.0 * c[2:s-1, s:s])) / (1.0+r) # Y:s
    f[1:1, 2:s-1] = (2.0 * c[1:1, 2:s-1] + (r-1.0) * p[1:1, 2:s-1] + 2.0*r*r*(c[2:2, 2:s-1] + c[1:1, 3:s] + c[1:1, 1:s-2] - 3.0 * c[1:1, 2:s-1])) / (1.0+r) # X:1
    f[s:s, 2:s-1] = (2.0 * c[s:s, 2:s-1] + (r-1.0) * p[s:s, 2:s-1] + 2.0*r*r*(c[s-1:s-1, 2:s-1] + c[s:s, 3:s] + c[s:s, 1:s-2] - 3.0 * c[s:s, 2:s-1])) / (1.0+r) # Y:s

    return 0
end

function prime_wave2d()
    speed = 10         # Propagation speed
    s = 16             # Array size (spatial resolution of the simulation)

    s2 = div(s, 2)
    s4 = div(s, 4)
    s8 = div(s, 8)
    s16 = div(s, 16)

    p = zeros(s, s) # past
    c = zeros(s, s) # current
    f = zeros(s, s) # future

    dt = 0.0001     # Time resolution of the simulation
    dx = 0.01       # Distance between elements
    r = speed * dt / dx

    n = 300

    for i = s2 - s16 : s2 + s16
        # Initial conditions
        c[i, s2 - s16 : s2 + s16] = - 2 * cos.(0.5 * 2 * pi / (s8) * (s2 - s16 : s2 + s16)) * cos.(0.5 * 2 * pi / (s8) * i)
        p[i, 1:s] = c[i, 1:s]
    end

    runWaveStencil(p, c, f, r, 0.0, n, s2, s4, s)
end

function wave2d(demo::Bool)

    speed = 10         # Propagation speed
    s = 512            # Array size (spatial resolution of the simulation)

    if demo
        stopTime = 0.1 # Time step at which to stop the main loop
    else
        stopTime = 0.05
    end

    s2 = div(s, 2)
    s4 = div(s, 4)
    s8 = div(s, 8)
    s16 = div(s, 16)

    p = zeros(s, s) # past
    c = zeros(s, s) # current
    f = zeros(s, s) # future

    dt = 0.0001     # Time resolution of the simulation
    dx = 0.01       # Distance between elements
    r = speed * dt / dx

    n = 300

    for i = s2 - s16 : s2 + s16
        # Initial conditions
        c[i, s2 - s16 : s2 + s16] = - 2 * cos.(0.5 * 2 * pi / (s8) * (s2 - s16 : s2 + s16)) * cos.(0.5 * 2 * pi / (s8) * i)
        p[i, 1:s] = c[i, 1:s]
    end

    # Main loop
    for t = 0 : dt : stopTime
        # Wave equation
        runWaveStencil(p, c, f, r, t, n, s2, s4, s)

        # Transparent corner handling
        f[1:1, 1:1] = (2 * c[1:1, 1:1] + (r-1) * p[1:1, 1:1] + 2*r*r* (c[2:2, 1:1] + c[1:1, 2:2] - 2*c[1:1, 1:1])) / (1+r) # X:1; Y:1
        f[s:s, s:s] = (2 * c[s:s, s:s] + (r-1) * p[s:s, s:s] + 2*r*r* (c[s-1:s-1, s:s] + c[s:s, s-1:s-1] - 2*c[s, s])) / (1+r) # X:s; Y:s
        f[1:1, s:s] = (2 * c[1:1, s:s] + (r-1) * p[1:1, s:s] + 2*r*r* (c[2:2, s:s] + c[1:1, s-1:s-1] - 2*c[1:1, s:s])) / (1+r) # X:1; Y:s
        f[s:s, 1:1] = (2 * c[s:s, 1:1] + (r-1) * p[s:s, 1:1] + 2*r*r* (c[s-1:s-1, 1:1] + c[s:s, 2:2] - 2*c[s:s, 1:1])) / (1+r) # X:s; Y:1

        # Rotate buffers for next iteration
        tmp = p
        p = c
        c = f
        f = tmp

        if demo
            if mod(t/dt, 10) == 0
                plot = Winston.imagesc(c)
                Winston.display(plot)
            end
        end
    end
end

function main()
    doc = """wave-2d.jl

Two-dimensional wave equation solver.

Usage:
  wave-2d.jl -h | --help
  wave-2d.jl [--demo]

Options:
  -h --help  Show this screen.
  --demo     Use settings that look good for the purposes of a demo, and plot the results using Winston.
"""
    arguments = docopt(doc)
    demo = arguments["--demo"]

    tic()
    prime_wave2d()
    println("SELFPRIMED ", toq())

    tic()
    wave2d(demo)
    println("SELFTIMED ", toq())
end

main()


text:
using Images
using ImageIO
using ImageMagick
using LinearAlgebra
using Pipe: @pipe
using ProgressMeter
using CSV
using DataFrames
X = load("jpgtest.jpg")

a = @pipe X |> Float64.(Gray.(_)) |> svd(_).S[1:10]'

data = Float64.(Gray.(X))
abssvd = svd(data)abs
a = Array{Float64}(undef, 0, 10)
@time a = vcat(a,re)

"""
Usage:
```julia
    recodeimage(path_to_folder_with_images, no_of_required_singular_vals)
```
Currently only supported for jpeg/png files
"""
function recodeimage(pathtoimage, n_singularvlas)
    imagelist = glob("*.jpg", pathtoimage)
    append!(imagelist,  glob("*.png", pathtoimage))
    recodedArray = Array{Float64}(undef, 0, n_singularvlas)
    @showprogress for images in imagelist
        X = load(images)
        img_singluar = @pipe X |> Float64.(Gray.(_)) |> svd(_).S[1:n_singularvlas]'
        recodedArray = vcat(recodedArray, img_singluar)
    end
    filename = joinpath(pathtoimage, "image_recoded.csv")
    CSV.write(filename,  DataFrame(recodedArray), writeheader=true)
end
using StatsBase
recodeimage(pwd(),2)
using WAV
using Pipe: @pipe
using MFCC
using Glob
using WAV

temp = Array{Float64}(undef,0,10)
x,fs = wavread("ENG_M.wav")
a = collect(mfcc(x,fs, numcep = 10)[1])
temp = vcat(temp,mean(a,dims=1))
using CSV
using DataFrames
using ProgressMeter
function recodeaudio(filepath,ncep)
    audiolist = glob("*.wav",filepath)
    recodedArray = Array{Float64}(undef, 0, ncep)
    @showprogress for files in audiolist
        x,fs = wavread(files)
        cep = collect(mfcc(x,fs, numcep = ncep)[1])# this is a matrix with ncep columns
        recodedArray = vcat(recodedArray, mean(cep,dims=1))
    end
    filename = joinpath(filepath, "audio_recoded.csv")
    CSV.write(filename,  DataFrame(recodedArray), writeheader=true)
end

using TextAnalysis
using Taro
using Languages
Taro.init()
function recodetext(pathtotxt, n_singularvals)
    files = glob("*.pdf", pathtotxt)
    docs = Array{StringDocument{String},1}[]
    getTitle(t) = TextAnalysis.sentence_tokenize(Languages.English(), t)[1]
   @showprogress for i in 1:length(files)
        meta,txt = Taro.extract(files[i]);
        txt = replace(txt, '\n' => ' ')
        title = getTitle(txt)
        dm = TextAnalysis.DocumentMetadata(Languages.English(), title, "", meta["Creation-Date"] )
        doc = StringDocument(txt, dm)
        docs = vcat(docs, doc)
    end
    crps = Corpus(docs)

    prepare!(crps, strip_non_letters | strip_punctuation | strip_case | strip_stopwords)
    stem!(crps)
    update_lexicon!(crps)
    update_inverse_index!(crps)
    U,S,Vt = @pipe crps |> DocumentTermMatrix(_) |> dtm(_, :dense) |> svd(_)
    if n_singularvals > length(S)
        println("You have set more number of singluar values than present \n so seting number of singluar values = length(singular values)")
        n_singularvals = length(S)
    end
    DTM = U[:,1:n_singularvals]*diagm(S[1:n_singularvals])*Vt[:,1:n_singularvals]'
    filename = joinpath(pathtotxt, "reconstructed_document_term_matrix.csv")
    CSV.write(filename,  DataFrame(DTM), writeheader=true)
end

a = glob("*.pdf", pwd())
docs= Any[]
getTitle(t) = TextAnalysis.sentence_tokenize(Languages.English(), t)[1]
for i in 1:length(a)
        meta,txt = Taro.extract(a)
        txt = replace(txt, '\n' => ' ')
        title = getTitle(txt)
        dm = TextAnalysis.DocumentMetadata(Languages.English(), title, "", meta["Creation-Date"] )
        doc = StringDocument(txt, dm)
        push!(docs, doc)

    end
end
a = recodetext(pwd(),50)

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
"""
# Summary SwitchPath
    mutable struct SwitchPath <: Scaler
# Fields
    slope::AbstractFloat
    scale::VarOrNil
Applies x -> x/(1+exp(-slope*scale) over a N-dimensional input x,
where slope is user defined, scale is a scalar Variable.
"""
mutable struct SwitchPath <: Scaler
    slope::AbstractFloat
    scale::VarOrNil
    function SwitchPath(slope::AbstractFloat; ndims::Int, type::Type=Array{Float32})
        @assert ndims >= 1 "ndims >= 1 shall be met, but got ndims=$ndims"
        shape = ntuple(i->1, ndims)
        scale = Variable{type}(Zeros(type, shape) .+ eltype(type)(7/slope), true, true, true);
        new(slope, scale)
    end
    function SwitchPath()
        new(1.0f0, nothing)
    end
end

function clone(this::SwitchPath; type::Type=Array{Float32})
    cloned = SwitchPath()
    cloned.slope = this.slope
    cloned.scale = clone(this.scale, type=type)
    return cloned
end

function Base.show(io::IO, m::SwitchPath)
    print(io, "SwitchPath(slope=$slope, scale=$(m.scale.value[1]); type=$TYPE)")
end

function paramsof(m::SwitchPath)
    params = Vector{Variable}(undef,1)
    params[1] = m.scale
    return params
end

function xparamsof(m::SwitchPath)
    xparams = Vector{XVariable}(undef,1)
    xparams[1] = ('w', m.scale)
    return xparams
end

function nparamsof(m::SwitchPath)
    return 1
end

function bytesof(m::SwitchPath, unit::String="MB")
    return blocksize(sizeof(m.scale), uppercase(unit))
end

function forward(m::SwitchPath, x::Variable{T}) where T
    k = m.slope  # f(x) = σ(k*x)
    a = m.scale  # scale param
    G = 1 / (1 + exp(-k * ᵛ(a)))
    y = Variable{T}(ᵛ(x) .* G, x.backprop)

    if x.backprop
        function SwitchPathBackward()
            if need2computeδ!(x) δ(x) .+=     δ(y) .* G                      end
            if need2computeδ!(a) δ(a) .+= sum(δ(y) .* k .* (1 .- G) .* ᵛ(y)) end
            ifNotKeepδThenFreeδ!(y);
        end
        push!(graph.backward, SwitchPathBackward)
    end
    return y
end

function predict(m::SwitchPath, x::AbstractArray)
    k = m.slope
    a = m.scale
    C = 1 / (1 + exp(-k * ᵛ(a)))
    return ᵛ(x) .* C
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
@RigidBodyDynamics.indextype BristleID

abstract type FrictionModel end

struct Bristle <: FrictionModel
    BristleID::BristleID
    τ::Float64
    k̄::Float64
	μs::Float64
	μd::Float64
	T̄s_μs::Float64
	T̄s_μd::Float64
	magic::Float64
	function Bristle(bristle_ID::BristleID; τ::Float64, k̄::Float64, μs::Float64, μd::Float64, magic::Float64)
		μs, μd = determine_μs_μd(μs, μd)
		T̄s_μs = 2 * μs
		T̄s_μd = 3 * μs
		return new(bristle_ID, τ, k̄, μs, μd, T̄s_μs, T̄s_μd, magic)
	end
end

struct Regularized <: FrictionModel
	μs::Float64
	μd::Float64
	v_c::Float64
	v_μs::Float64
	v_μd::Float64
    function Regularized(v_c; μs::Float64, μd::Float64)
		μs, μd = determine_μs_μd(μs, μd)
		v_μs = 2 * v_c
		v_μd = 3 * v_c
		return new(μs, μd, v_c, v_μs, v_μd)
	end
end

struct ContactInstructions
    id_1::MeshID
    id_2::MeshID
    χ::Float64
    FrictionModel::Union{Regularized,Bristle}
	quad::TriTetQuadRule{3,N} where {N}
    function ContactInstructions(id_tri::MeshID, id_tet::MeshID, fric_model::Union{Regularized,Bristle}; χ::Float64,
			n_quad_rule::Int64)

		(1 <= n_quad_rule <= 2) || error("only quadrature rules 1 (first order) and 2 (second? order) are currently implemented")
		quad = getTriQuadRule(n_quad_rule)
        return new(id_tri, id_tet, χ, fric_model, quad)
    end
end

struct TractionCache{T}
	n̂::SVector{3,T}
    r_cart::SVector{3,T}
    dA::T
    p::T
end

@inline calc_p_dA(t::TractionCache) = t.p * t.dA

mutable struct spatialStiffness{T}
	K::Hermitian{T,Matrix{T}}
	K̄::Hermitian{T,Matrix{T}}
    σ_sqrt::Diagonal{T,MVector{6,T}}
    K̄⁻¹_sqrt::MMatrix{6,6,T,36}
	S⁻¹::Diagonal{T,MVector{6,T}}
	mul_pre::MMatrix{6,6,T,36}
    function spatialStiffness{T}() where {T}
		K = Hermitian(Matrix(T.(zeros(6,6) .+ NaN)))
		K̄ = Hermitian(Matrix(T.(zeros(6,6) .+ NaN)))
		σ_sqrt = Diagonal(MVector{6,T}(NaN, NaN, NaN, NaN, NaN, NaN))
        K̄⁻¹_sqrt = MMatrix{6,6,T,36}(T.(zeros(6,6) .+ NaN))
		S⁻¹ = Diagonal(MVector{6,T}(NaN, NaN, NaN, NaN, NaN, NaN))
        mul_pre = MMatrix{6,6,T,36}(T.(zeros(6,6)) .+ NaN)
        return new(K, K̄, σ_sqrt, K̄⁻¹_sqrt, S⁻¹, mul_pre)
    end
end

mutable struct TypedElasticBodyBodyCache{T}
    spatialStiffness::spatialStiffness{T}
	TractionCache::VectorCache{TractionCache{T}}
    quad::TriTetQuadRule{3,N} where {N}
    mesh_1::MeshCache
    mesh_2::MeshCache
    x_rʷ_r²::Transform3D{T}
    x_r²_rʷ::Transform3D{T}
    x_r¹_r²::Transform3D{T}
	x_r²_r¹::Transform3D{T}
    twist_r²_r¹_r²::Twist{T}
    χ::Float64
    Ē::Float64
    wrench::Wrench{T}
    function TypedElasticBodyBodyCache{T}() where {T}
        trac_cache = VectorCache{TractionCache{T}}()
        K = spatialStiffness{T}()
        return new{T}(K, trac_cache)
    end
end

struct TypedMechanismScenario{T}
    state::MechanismState{T}
    s::SegmentedVector{BristleID,T,Base.OneTo{BristleID},Array{T,1}}
    result::DynamicsResult{T}
    ṡ::SegmentedVector{BristleID,T,Base.OneTo{BristleID},Array{T,1}}
    f_generalized::Vector{T}
    bodyBodyCache::TypedElasticBodyBodyCache{T}
    GeometricJacobian::RigidBodyDynamics.CustomCollections.CacheIndexDict{BodyID,Base.OneTo{BodyID},Union{Nothing,GeometricJacobian{Array{T,2}}}}
    torque_third_law::Vector{T}
    τ_ext::Vector{T}
    rhs::Vector{T}
    nv::Int64

    function TypedMechanismScenario{T}(mechanism::Mechanism, v_path, body_ids, n_bristle_pairs::Int64) where {T}

        function makeJacobian(v_path, state::MechanismState{T}, body_ids::Base.OneTo{BodyID}) where {T}
            v_jac = RigidBodyDynamics.BodyCacheDict{Union{Nothing,GeometricJacobian{Array{T,2}}}}(body_ids)
            fill_with_nothing!(v_jac)
            for (body_id_k, path_k) = v_path
                if path_k != nothing  # body has no meshes attached to it
                    if body_id_k != BodyID(root_body(state.mechanism))
                        new_jac = geometric_jacobian(state, path_k)
                        fill!(new_jac.linear, NaN + zero(T))
                        fill!(new_jac.angular, NaN + zero(T))
                        v_jac[body_id_k] = new_jac
                    end
                end
            end
            return v_jac
        end
        function_Int64_six(k) = 6

        state = MechanismState{T}(mechanism)
        result = DynamicsResult{T}(mechanism)
        f_generalized = Vector{T}(undef, num_positions(mechanism))
        bodyBodyCache = TypedElasticBodyBodyCache{T}()
        v_jac = makeJacobian(v_path, state, body_ids)
        n_dof_bristle = 6 * n_bristle_pairs
        s = SegmentedVector{BristleID}(zeros(T, n_dof_bristle), Base.OneTo(BristleID(n_bristle_pairs)), function_Int64_six)
        ṡ = SegmentedVector{BristleID}(zeros(T, n_dof_bristle), Base.OneTo(BristleID(n_bristle_pairs)), function_Int64_six)
        torque_third_law = Vector{T}(undef, num_positions(mechanism))
        τ_ext = zeros(T, num_positions(mechanism))
        rhs = zeros(T, num_positions(mechanism))
        nv = num_velocities(mechanism)
        return new{T}(state, s, result, ṡ, f_generalized, bodyBodyCache, v_jac, torque_third_law, τ_ext, rhs, nv)
    end

	function TypedMechanismScenario{T}(mechanism::Mechanism) where {T}
		return new{T}(MechanismState{T}(mechanism))
	end
end

mutable struct DiscreteControl
	dt::Float64
	t_last::Float64
	control::Function
	DiscreteControl(control::Function, dt=0.01) = new(dt, 0.0, control)
end

"""
$(TYPEDEF)

A `MechanismScenario` contains contact information for an entire mechanism.

Type parameters:
* `T`: the `ForwardDiff.Dual` scalar type
"""
mutable struct MechanismScenario{T}
	float::TypedMechanismScenario{Float64}
	dual::TypedMechanismScenario{T}
    body_ids::Base.OneTo{BodyID}
    mesh_ids::Base.OneTo{MeshID}
    bristle_ids::Base.OneTo{BristleID}
	MeshCache::RigidBodyDynamics.CustomCollections.CacheIndexDict{MeshID,Base.OneTo{MeshID},MeshCache}
	ContactInstructions::Vector{ContactInstructions}
	TT_Cache::TT_Cache
	de::Function
	continuous_controller::Union{Nothing,Function}
	discrete_controller::Union{Nothing,DiscreteControl}
    τ_ext::Vector{Float64}
    path::RigidBodyDynamics.CustomCollections.IndexDict{BodyID,Base.OneTo{BodyID},Union{Nothing,RigidBodyDynamics.Graphs.TreePath{RigidBody{Float64},Joint{Float64,JT} where JT<:JointType{Float64}}}}

	function MechanismScenario(; de::Function=calcXd!, N_chunk::Int64=6,
			continuous_controller::Union{Nothing,Function}=nothing,
			discrete_controller::Union{Nothing,DiscreteControl}=nothing,
			gravity::SVector{3,Float64}=SVector{3,Float64}(0.0, 0.0, -round(9.8054, digits=8, base=2) ))
		#
		mechanism = Mechanism(RigidBody{Float64}("world"); gravity=gravity)  # create empty mechanism
        T = Dual{Nothing,Float64,N_chunk}
		cache_float = TypedMechanismScenario{Float64}(mechanism)
		cache_dual = TypedMechanismScenario{T}(mechanism)
		body_ids = Base.OneTo(BodyID(1))
		mesh_ids = Base.OneTo(MeshID(0))
		bristle_ids = Base.OneTo(BristleID(0))
		mesh_cache = MeshCacheDict{MeshCache}(mesh_ids)
		c_ins = Vector{ContactInstructions}()

		return new{T}(cache_float, cache_dual, body_ids, mesh_ids, bristle_ids, mesh_cache, c_ins, TT_Cache(), de,
			continuous_controller, discrete_controller)
    end
end

"""
$(SIGNATURES)

Calculates `MechanismScenario` information needed for simulation.
"""
function finalize!(m::MechanismScenario{T}) where {T}
	function makePaths(mechanism::Mechanism, mesh_cache::MeshCacheDict{MeshCache}, body_ids::Base.OneTo{BodyID})
		the_type = Union{Nothing,RigidBodyDynamics.Graphs.TreePath{RigidBody{Float64},Joint{Float64,JT} where JT<:JointType{Float64}}}
		v_path = RigidBodyDynamics.BodyDict{the_type}(body_ids)
		fill_with_nothing!(v_path)
		for mesh_k = values(mesh_cache)
			body_id = mesh_k.BodyID
			if body_id != BodyID(root_body(mechanism))
				v_path[body_id] = path(mechanism, root_body(mechanism), bodies(mechanism)[body_id])
			end
		end
		return v_path
	end

	# quad = getTriQuadRule(ifelse(NQ == 1, 1, 2))
	mechanism = get_mechanism(m)
	body_ids = Base.OneTo(last(bodies(mechanism)).id)
	m.body_ids = body_ids
	n_bristle_pairs = length(m.bristle_ids)
	m.τ_ext = zeros(Float64, num_positions(mechanism))
	m.path = makePaths(mechanism, m.MeshCache, body_ids)
	cache_path = m.path
	m.float = TypedMechanismScenario{Float64}(mechanism, cache_path, body_ids, n_bristle_pairs)
	m.dual  = TypedMechanismScenario{T}(mechanism, cache_path, body_ids, n_bristle_pairs)
	return nothing
end

num_partials(m::MechanismScenario{Dual{Nothing,Float64,N_partials}}) where {N_partials} = N_partials
function num_x(m::MechanismScenario{T}) where {T}
    tm = m.float
    mechanism = tm.state.mechanism
    return num_positions(mechanism) + num_velocities(mechanism) + length(tm.s)
 end
type_dual(m::MechanismScenario{T}) where {T} = T

function get_state(m::MechanismScenario)
    x = zeros(num_x(m))
    copyto!(x, m.float)
    return x
end

function set_state_spq!(m::MechanismScenario, joint::Joint; rot::Rotation=one(Quat{Float64}),
        trans::SVector{3,Float64}=zeros(SVector{3,Float64}), w::SVector{3,Float64}=zeros(SVector{3,Float64}),
        vel::SVector{3,Float64}=zeros(SVector{3,Float64}))

    rot = components(SPQuat(rot))
    state = m.float.state
    set_configuration!(state, joint, vcat(rot, trans))
    set_velocity!(state, joint, vcat(w, vel))
    return nothing
end

function addMesh!(m::MechanismScenario, mesh::MeshCache)
    mesh_ids_old = m.mesh_ids
    mesh_ids_new = Base.OneTo(MeshID(length(m.mesh_ids) + 1))
    mesh_cache = MeshCacheDict{MeshCache}(mesh_ids_new)
    for id = mesh_ids_old
        mesh_cache[id] = m.MeshCache[id]
    end
    mesh_cache[mesh_ids_new[end]] = mesh
    m.MeshCache = mesh_cache
    m.mesh_ids = mesh_ids_new
    return nothing
end

get_mechanism(m::MechanismScenario) = m.float.state.mechanism

"""
$(SIGNATURES)

Creates a body with the contact geometry and inertia specified by the inputs.
Internally, this function calls [`add_contact!`](@ref) and [`add_body!`](@ref).
"""
function add_body_contact!(m::MechanismScenario, name::String, e_mesh::eMesh;
        i_prop::InertiaProperties,
        c_prop::Union{Nothing,ContactProperties}=nothing,
        body::Union{RigidBody{Float64},Nothing}=nothing,
        joint::JT=SPQuatFloating{Float64}(),
        dh::basic_dh=one(basic_dh{Float64})) where {JT<:JointType}

    nt = add_body!(m, name, e_mesh, i_prop=i_prop, body=body, joint=joint, dh=dh)
    nt_new_contact = add_contact!(m, name, e_mesh, c_prop=c_prop, body=nt.body)
    return NamedTuple{(:body, :joint, :id)}((nt.body, nt.joint, nt_new_contact.id))
end

"""
$(SIGNATURES)

Adds contact geometry without creating a new body.
In addition to being used by other functions, this function allows you to add contact geometry to an imported mechanism,
or add multiple contact geometries to a single body.
"""
function add_contact!(m::MechanismScenario, name::String, e_mesh::eMesh;  c_prop::Union{Nothing,ContactProperties}=nothing,
        body::Union{RigidBody{Float64},Nothing}=nothing)

	function verify_eMesh_ContactProperties(eM::eMesh{T1,T2}, c_prop::Union{Nothing,ContactProperties}) where {T1,T2}
		(T1 == Tri) && (T2 == Tet)         && error("eMesh has triangles and tets. Use as_tri_eMesh or as_tet_eMesh to convert eMesh.")
		(T1 == nothing) && (T2 == nothing) && error("eMesh has neither triangles and tets.")
		(T1 == Tri) && (c_prop != nothing) && error("Using ContactProperties for triangular eMesh")
		(T2 == Tet) && (c_prop == nothing) && error("Using nothing as ContactProperties for tet eMesh")
	end

    verify_eMesh_ContactProperties(e_mesh, c_prop)
    body = return_body_never_nothing(get_mechanism(m), body)
    tree = eMesh_to_tree(e_mesh)
	mesh = MeshCache(name, e_mesh, tree, body, c_prop)
    addMesh!(m, mesh)
    return NamedTuple{(:id,)}((find_mesh_id(m, mesh),))
end

"""
$(SIGNATURES)

Adds a body to the mechanism with whose inertia properties are those of the contact geometry in `e_mesh` and inertia
properties in `i_prop`.
In addition to being used by other functions, this function is helpful to create bodies with specific inertia properties
in an intuitive manner.
"""
function add_body!(m::MechanismScenario, name::String, e_mesh::eMesh; i_prop::InertiaProperties,
        body::Union{RigidBody{Float64},Nothing}=nothing, joint::JT=SPQuatFloating{Float64}(),
        dh::basic_dh=one(basic_dh{Float64})) where {JT<:JointType}

    mesh_inertia_info = makeInertiaInfo(e_mesh, i_prop)
    return add_body_from_inertia!(get_mechanism(m), name, mesh_inertia_info, joint=joint, body=body, dh=dh)
end

return_body_never_nothing(mechanism::Mechanism, body::Nothing) = root_body(mechanism)
return_body_never_nothing(mechanism::Mechanism, body::RigidBody{Float64}) = body

function add_body_from_inertia!(mechanism::Mechanism, name::String, mesh_inertia_info::MeshInertiaInfo;
        joint::JT=SPQuatFloating{Float64}(), body::Union{RigidBody{Float64},Nothing}=nothing,
        dh::basic_dh{Float64}=one(basic_dh{Float64})) where {JT<:JointType}
    # TODO: check that a spherical floating joint isn't added

    body_parent = return_body_never_nothing(mechanism, body)
    body_child = newBodyFromInertia(name, mesh_inertia_info)
    j_parent_child, x_parent_child = outputJointTransform_ParentChild(body_parent, body_child, joint, dh)
    attach!(mechanism, body_parent, body_child, j_parent_child, joint_pose=x_parent_child)
    return NamedTuple{(:body, :joint)}((body_child, j_parent_child))
end

default_χ() = 0.5
default_μ() = 0.3

determine_μs_μd(μs::Nothing, μd::Nothing) = (default_χ(), default_χ())
determine_μs_μd(us::Float64, μd::Nothing) = error("need to specify μd")
determine_μs_μd(μs::Nothing, μd::Float64) = (μd, μd)
function determine_μs_μd(μs::Float64, μd::Float64)
	(μd <= μs) || error("something is wrong")
	return μs, μd
end

# TODO: explain what the tunable paramaters are.
"""
$(SIGNATURES)

Adds regularized friction to a scenario.
Returns the added contact instruction.
"""
function add_friction_regularize!(m::MechanismScenario, mesh_id_1::MeshID, mesh_id_2::MeshID;
		μs::Union{Float64,Nothing}=nothing,
		μd::Union{Float64,Nothing}=nothing,
		χ::Float64=default_χ(), v_tol::Float64=0.01,
		n_quad_rule::Int64=2
		)
	#
	μs, μd = determine_μs_μd(μs, μd)
    regularized = Regularized(v_tol, μs=μs, μd=μd)
    return add_friction!(m, mesh_id_1, mesh_id_2, regularized, χ=χ, n_quad_rule=n_quad_rule)
end

# TODO: explain what the tunable paramaters are.
"""
$(SIGNATURES)

Adds bristle friction to a scenario.
Returns the added contact instruction.
"""
function add_friction_bristle!(m::MechanismScenario, mesh_id_1::MeshID, mesh_id_c::MeshID; τ::Float64=0.05, k̄=1.0e4,
		μs::Union{Float64,Nothing}=nothing,
		μd::Union{Float64,Nothing}=nothing,
		χ::Float64=default_χ(),
		n_quad_rule::Int64=2,
		magic::Float64=1.0e-3)
	#
	μs, μd = determine_μs_μd(μs, μd)
    isa(μd, Nothing) || (0 < μd) || error("μd cannot be 0 for bristle friction")
    bristle_id = BristleID(1 + length(m.bristle_ids))
    bf = Bristle(bristle_id, τ=τ, k̄=k̄, μs=μs, μd=μd, magic=magic)
    m.bristle_ids = Base.OneTo(bristle_id)
    return add_friction!(m, mesh_id_1, mesh_id_c, bf, χ=χ, n_quad_rule=n_quad_rule)
end

add_friction!(m::MechanismScenario, id_1::MeshID, id_2::MeshID, fric_model::FrictionModel; χ::Float64,
	n_quad_rule::Int64) = add_friction!(m, id_1, id_2, m.MeshCache[id_1], m.MeshCache[id_2], fric_model, χ=χ, n_quad_rule=n_quad_rule)

function add_friction!(m::MechanismScenario, id_1::MeshID, id_2::MeshID, m_1::MeshCache{Nothing,Tet},
		m_2::MeshCache{Tri,Nothing}, fric_model::Union{Regularized,Bristle}; χ::Float64,
		n_quad_rule::Int64)

	return add_friction!(m, id_2, id_1, fric_model; χ=χ, n_quad_rule=n_quad_rule)
end

function add_friction!(m::MechanismScenario, id_1::MeshID, id_2::MeshID, m_1::MeshCache{T1,T2},
		m_2::MeshCache{Nothing,Tet}, fric_model::Union{Regularized,Bristle}; χ::Float64,
		n_quad_rule::Int64) where {T1,T2}

	c_ins_new = ContactInstructions(id_1, id_2, fric_model, χ=χ, n_quad_rule=n_quad_rule)
	push!(m.ContactInstructions, c_ins_new)
	return c_ins_new
end


text:
using Documenter, AbstractGPs

makedocs(;
    modules=[AbstractGPs],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/JuliaGaussianProcesses/AbstractGPs.jl/blob/{commit}{path}#L{line}",
    sitename="AbstractGPs.jl",
    authors="willtebbutt <wt0881@my.bristol.ac.uk>",
    assets=String[],
)

deploydocs(;
    repo="github.com/JuliaGaussianProcesses/AbstractGPs.jl",
)


text:
using ExperimentalSetup
const ES = ExperimentalSetup
using Base.Test, DataStructures
import Base: ==
==(a::Base.KeyIterator, b::Base.KeyIterator) = length(a)==length(b) && all(k->in(k,b), a)
==(a::ES.Log, b::ES.Log) = a.last == b.last && a.md == b.md && a.reps == b.reps
==(a::ES.Metadata, b::ES.Metadata) = a.factors == b.factors && a.nfactors == b.nfactors && a.setups == b.setups
==(a::ES.Factor, b::ES.Factor) = a.name == b.name && a.levels == b.levels
==(a::Vector{ES.Factor}, b::Vector{ES.Factor}) = length(a) == length(b) && all(a[i] == b[i] for i in 1:length(a))
==(a::ES.Setup, b::ES.Setup) = a.levels == b.levels && a.n.x == b.n.x
==(a::Vector{ES.Setup}, b::Vector{ES.Setup}) = length(a) == length(b) && all(a[i] == b[i] for i in 1:length(a))
==(a::ES.Rep, b::ES.Rep) = a.comment == b.comment && a.replicate == b.replicate && a.setup == b.setup
==(a::Vector{ES.Rep}, b::Vector{ES.Rep}) = length(a) == length(b) && all(a[i] == b[i] for i in 1:length(a))
function create_log(f=string.(1:5))
    factors = [Factor(f[i], f[1:i]) for i in linearindices(f)]
    a = ES.Log(ES.Metadata(factors))
    push!(a, repeat(f[1,:], outer=[length(f)]), "comment")
    push!(a, f, "comment")
end

@testset "other" begin

    a = create_log()
    @test a.md.nfactors == 5
    @test 1 ∈ a

    empty!(a)
    @test isempty(a.reps)
    @test isempty(a.md.setups)

    levels = [string(rand(1:i)) for i in 1:5]
    comment = "some comment"
    push!(a, levels, comment)
    @test ES.pop(a, a.last) == (levels, comment)

end

@testset "push" begin

    b = create_log()
    # replicate
    push!(b, repeat(["1"], outer=[5]), "comment")
    @test b.reps[end].setup == 1
    @test b.reps[end].replicate == 2

    # comment
    push!(b, repeat(["1"], outer=[5]), "different comment")
    @test b.reps[end].setup == 1
    @test b.reps[end].replicate == 3

    # new
    push!(b, [i > 1 ? "2" : "1" for i in 1:5], "comment")
    @test b.reps[end].setup == 3
    @test b.reps[end].replicate == 1

    # errors
    @test_throws AssertionError push!(b, repeat(["1"], outer=[6]), "comment")
    @test_throws AssertionError push!(b, repeat(["2"], outer=[5]), "comment")

end

@testset "delete" begin

    # delete a member
    b = create_log()
    delete!(b, 1)
    @test 1 ∉ b

    # delete something that isn't there
    @test_throws KeyError delete!(b, 1)

    # delete a factor
    b = create_log()
    c = delete(b, b.md.factors[2])
    @test b.md.factors[2] ∉ c.md.factors
    @test all(length(f.levels) == 4 for f in c.md.setups)
    @test all(f.levels[2] ≠ "2" for f in c.md.setups)

    # delete a level
    b = create_log()
    push!(b, [i > 1 ? "2" : "1" for i in 1:5], "comment")
    c = delete(b, b.md.factors[5], "2")
    f = 1:5
    factors = [i ≠ 5 ? Factor("$i", ["$j" for j in 1:i]) : Factor("$i", ["$j" for j in filter(x -> x ≠ 2, 1:i)]) for i in f]
    a = ES.Log(ES.Metadata(factors))
    push!(a, repeat(["1"], outer=[length(f)]), "comment")
    push!(a, ["$i" for i in f], "comment")
    push!(a, ["1", "2", "2", "2", "1"], "comment")

    @test c == a

end

@testset "replace" begin

    a = create_log()
    new = ([i < 3 ? "$i" : "3" for i in 1:5], "comment")
    old = 1
    replace!(a, old, new...)
    @test ES.pop(a, old) == new

end

@testset "combine" begin

    a = create_log()
    b = create_log(string.('a':'g'))
    c = ES.combine(a,b)
    @test c == ExperimentalSetup.Log(ExperimentalSetup.Metadata(ExperimentalSetup.Factor[ExperimentalSetup.Factor("1", String["1"]), ExperimentalSetup.Factor("2", String["1", "2"]), ExperimentalSetup.Factor("3", String["1", "2", "3"]), ExperimentalSetup.Factor("4", String["1", "2", "3", "4"]), ExperimentalSetup.Factor("5", String["1", "2", "3", "4", "5"]), ExperimentalSetup.Factor("a", String["a"]), ExperimentalSetup.Factor("b", String["a", "b"]), ExperimentalSetup.Factor("c", String["a", "b", "c"]), ExperimentalSetup.Factor("d", String["a", "b", "c", "d"]), ExperimentalSetup.Factor("e", String["a", "b", "c", "d", "e"]), ExperimentalSetup.Factor("f", String["a", "b", "c", "d", "e", "f"]), ExperimentalSetup.Factor("g", String["a", "b", "c", "d", "e", "f", "g"])], ExperimentalSetup.Setup[ExperimentalSetup.Setup([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], Base.RefValue{Int64}(2)), ExperimentalSetup.Setup([1, 2, 3, 4, 5, 1, 1, 1, 1, 1, 1, 1], Base.RefValue{Int64}(1)), ExperimentalSetup.Setup([1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7], Base.RefValue{Int64}(1))]), 0x0000000000000004, SortedDict(0x0000000000000001=>ExperimentalSetup.Rep(1, 1, "comment"),0x0000000000000002=>ExperimentalSetup.Rep(2, 1, "comment"),0x0000000000000003=>ExperimentalSetup.Rep(1, 2, "comment"),0x0000000000000004=>ExperimentalSetup.Rep(3, 1, "comment")))

end

@testset "json" begin
    const filename = "log.json"
    # create a log object to use for our tests
    log = create_log()

    # ExperimentalSetup.parseJSON tests: valid cases
    # if we convert a Log to JSON and back again, the resulting Log should be identical
    jsonString = JSON.json(log)
    @test log == ES.parseJSON(jsonString)
    # we should be able to parse an empty log without any errors
    @test ES.Log() == ES.parseJSON(JSON.json(ES.Log()))

    # ExperimentalSetup.parseJSON tests: invalid cases
    # we expect an exception if we try to parse an empty string
    @test_throws Exception ES.parseJSON("")
    # we also expect an exception if we try to parse an incomplete/malformed JSON string
    @test_throws Exception ES.parseJSON(jsonString[1:end-1])

    # ExperimentalSetup.saveJSON tests: valid cases
    # save the object to a file
    ES.saveJSON(filename, log)
    # we expect the file to have been created
    @test isfile(filename)
    open(filename, "r") do f
        # we expect the file to contain the same JSON
        @test JSON.json(log) == readstring(f)
    end
    # we expect an error if we pass a filename that exists with overwrite = false
    @test_throws Exception ES.saveJSON(filename, log, false)

    # ExperimentalSetup.saveJSON tests: invalid cases
    # we expect an error if we pass an empty filename
    @test_throws Exception ES.saveJSON("", log)
    # we expect an error if we pass an empty Log
    @test_throws Exception ES.saveJSON(filename, None)

    # ExperimentalSetup.loadJSON tests: valid cases
    # we expect the loaded object to be the same as the saved object
    parsedLogFromFile = ES.loadJSON(filename)
    @test log == parsedLogFromFile

    # ExperimentalSetup.loadJSON tests: invalid cases
    # we expect an error if we pass an empty filename
    @test_throws Exception ES.loadJSON("")
    # we expect an error if we pass a file that does not exist
    @test_throws Exception ES.loadJSON("this.file.should.not.exist")
    # we expect an error if we pass a file that does not contain valid JSON
    open(filename, "w") do f
        write(f, "This is not a JSON string!")
    end
    @test_throws Exception parsedLogFromFile = ES.loadJSON(filename)

    # cleanup
    rm(filename)
end


text:
"""
    Aquifer

Abstract type representing an aquifer, either confined or unconfined.

The vertically averaged governing equation of an unconfined, inhomogeneous and
isotropic aquifer in one dimension can be written as:

S * ∂ϕ / ∂t = ∂ / ∂x * (kH * ∂ϕ / ∂x) + Q

with:
* S: storativity (or storage coefficient)
* ϕ: hydraulic head (groundwater level)
* t: time
* k: conductivity
* H: H the (saturated) aquifer height: groundwater level - aquifer bottom
  elevation
* η: elevation of aquifer bottom
* Q: fluxes from boundary conditions (e.g. recharge or abstraction)

The simplest finite difference formulation is forward in time, central in space,
and can be written as:

Sᵢ * (ϕᵢᵗ⁺¹ - ϕᵢᵗ) / Δt = -Cᵢ₋₁ * (ϕᵢ₋₁ - ϕᵢ) -Cᵢ * (ϕᵢ₊₁ - ϕᵢ) + Qᵢ

with:
* ᵢ as cell index
* ᵗ as time index
* Δt as step size
* Cᵢ₋₁ as the intercell conductance between cell i-1 and i
* Cᵢ as the intercell conductance between cell i and i+1

Conductance is defined as:

C = kH * w / l

with:
* w the width of the cell to cell connection
* l the length of the cell to cell connection

k and H may both vary in space; intercell conductance is therefore an average
using the properties of two cells. See the documentation below.

There is only one unknown, ϕᵢᵗ⁺¹. Reshuffling terms:

ϕᵢᵗ⁺¹ = ϕᵢᵗ + (Cᵢ₋₁ * (ϕᵢ - ϕᵢ₋₁) + Cᵢ * (ϕᵢ₊₁ - ϕᵢ) + Qᵢ) * Δt / Sᵢ

This can be generalized to two dimensions, for both regular and irregular cell
connectivity.

See this paper for more details:
     Chu, W. S., & Willis, R. (1984). An explicit finite difference model for
     unconfined aquifers. Groundwater, 22(6), 728-734.

Boundary conditions can be classified into three categories:
* specified head (Dirichlet)
* specified flux (Neumann)
* head-dependent flux (Robin)

Neumann and Robin conditions are implemented by adding to or subtracting from a
net (lumped) cell flux. Dirichlet conditions are special cased, since they
cannot (easily) be implemented via the flux, but the head is set directly
instead.
"""
abstract type Aquifer end

abstract type AquiferBoundaryCondition end

"""
    ConfinedAquifer{T} <: Aquifer

Confined aquifers are overlain by a poorly permeable confining layer (e.g.
clay). No air can get in to fill the pore space so that the aquifer always
remains fully saturated. For a confined aquifer, water will always flow along
the complete height over the aquifer and transmissivity kH (m d⁻¹) is a
constant.

Specific storage is the amount of water an aquifer releases per unit change in
hydraulic head, per unit volume of aquifer, as the aquifer and the groundwater
itself is compressed. Its value is much smaller than specific yield, between
1E-5 (stiff) and 0.01 (weak).

NOTA BENE: **specific** storage is per m of aquifer (conf. specific weight).
**Storativity** or (**storage coefficient**) is for the entire aquifer (conf.
transmissivity).
"""
struct ConfinedAquifer{T} <: Aquifer
    head::Vector{T}  # hydraulic head [m]
    k::Vector{T}  # horizontal conductivity [m m⁻¹]
    top::Vector{T} # top of groundwater layer [m]
    bottom::Vector{T} # bottom of groundwater layer
    area::Vector{T} # area of cell
    specific_storage::Vector{T} # [m m⁻¹ m⁻¹]
    storativity::Vector{T} # [m m⁻¹]
    conductance::Vector{T} # Confined aquifer conductance is constant
end

"""
    UnconfinedAquifer{T} <: Aquifer

The upper boundary of an unconfined aquifer is the water table (the phreatic
surface).

Specific yield (or drainable porosity) represents the volumetric fraction the
aquifer will yield when all water drains and the pore volume is filled by air
instead. Specific yield will vary roughly between 0.05 (clay) and 0.45 (peat)
(Johnson, 1967).
"""
struct UnconfinedAquifer{T} <: Aquifer
    head::Vector{T}  # hydraulic head [m]
    k::Vector{T}  # horizontal conductivity [m d⁻¹]
    top::Vector{T} # top of groundwater layer [m]
    bottom::Vector{T} # bottom of groundwater layer
    area::Vector{T}
    specific_yield::Vector{T} # [m m⁻¹]
    conductance::Vector{T} #
    # Unconfined aquifer conductance is computed with degree of saturation
end

storativity(A::UnconfinedAquifer) = A.specific_yield
storativity(A::ConfinedAquifer) = A.storativity

"""
    harmonicmean_conductance(k1, k2, H1, H2, l1, l2, width)

The harmonic mean is the exact interblock transmissivity for steady-state
one-dimensional flow with no recharge if the transmissivity is assumed to be
spatially uniform over each finite difference block, changing abruptly at the
block interface.

Refer to:

    Goode, D. J., & Appel, C. A. (1992). Finite-Difference Interblock Transmissivity
    for Unconﬁned Aquifers and for Aquifers having Smoothly Varying Transmissivity.
    Water-resources investigations report, 92, 4124.
"""
function harmonicmean_conductance(k1, k2, H1, H2, l1, l2, width)
    kH1 = k1 * H1
    kH2 = k2 * H2
    if (kH1 * kH2) > 0.0
        return width * kH1 * kH2 / (kH1 * l2 + kH2 * l1)
    else
        return 0.0
    end
end

function saturated_thickness(aquifer::UnconfinedAquifer, index::Int)
    min(aquifer.top[index], aquifer.head[index]) - aquifer.bottom[index]
end

function saturated_thickness(aquifer::ConfinedAquifer, index::Int)
    aquifer.top[index] - aquifer.bottom[index]
end

"""
    horizontal_conductance(i, j, nzi, aquifer, C)

Compute fully saturated horizontal conductance for a single connection between two cells
(indexed with `i` and `j`). Geometry characteristics are taken from the
connectivity struct `C`, using the non-zero index (nzi) of its CSC data
structure.
"""
function horizontal_conductance(
    i::Int,
    j::Int,
    nzi::Int,
    aquifer::A,
    connectivity::Connectivity,
) where A <: Aquifer
    k1 = aquifer.k[i]
    k2 = aquifer.k[j]
    H1 = aquifer.top[i] - aquifer.bottom[i]
    H2 = aquifer.top[j] - aquifer.bottom[j]
    length1 = connectivity.length1[nzi]
    length2 = connectivity.length2[nzi]
    width = connectivity.width[nzi]
    return harmonicmean_conductance(k1, k2, H1, H2, length1, length2, width)
end

"""
    initialize_conductance!(aquifer::A, connectivity::Connectivity) where A <: Aquifer

Conductance for a confined aquifer is constant, and only has to be set once.
For an unconfined aquifer, conductance is computed per timestep by multiplying by
degree of saturation [0.0 - 1.0].
"""
function initialize_conductance!(aquifer::A, connectivity::Connectivity) where A <: Aquifer
    for i in 1:connectivity.ncell
        # Loop over connections for cell j
        for nzi in connections(connectivity, i)
            j = connectivity.rowval[nzi]
            aquifer.conductance[nzi] = horizontal_conductance(i, j, nzi, aquifer, connectivity)
        end
    end
end

function conductance(aquifer::ConfinedAquifer, i, j, nzi)
    return aquifer.conductance[nzi]
end

"""
    conductance(aquifer::UnconfinedAquifer, connectivity::Connectivity)

This computes the conductance for an unconfined aquifer using the "upstream
saturated fraction" as the MODFLOW documentation calls it. In this approach, the
saturated thickness of a cell-to-cell is approximated using the cell with the
highest head. This results in a consistent overestimation of the saturated
thickness, but it avoids complexities related with cell drying and rewetting,
such as having to define a "wetting threshold" or a "wetting factor".

See the documentation for MODFLOW-NWT or MODFLOW6 for more background:
Niswonger, R.G., Panday, Sorab, and Ibaraki, Motomu, 2011, MODFLOW-NWT, A Newton
formulation for MODFLOW-2005: U.S. Geological Survey Techniques and Methods
6-A37, 44 p.

Langevin, C.D., Hughes, J.D., Banta, E.R., Niswonger, R.G., Panday, Sorab, and
Provost, A.M., 2017, Documentation for the MODFLOW 6 Groundwater Flow Model:
U.S. Geological Survey Techniques and Methods, book 6, chap. A55, 197 p.,
https://doi.org/10.3133/tm6A55.

For background on drying and rewetting, see:
McDonald, M.G., Harbaugh, A.W., Orr, B.R., and Ackerman, D.J., 1991, A method of
converting no-flow cells to variable-head cells for the U.S. Geological Survey
modular finite-difference groundwater flow model: U.S. Geological Survey
Open-File Report 91-536, 99 p
"""
function conductance(aquifer::UnconfinedAquifer, i, j, nzi)
    ϕᵢ = aquifer.head[i]
    ϕⱼ = aquifer.head[j]
    if ϕᵢ >= ϕⱼ
        saturation = saturated_thickness(aquifer, i) / (aquifer.top[i] - aquifer.bottom[i])
    else
        saturation = saturated_thickness(aquifer, j) / (aquifer.top[j] - aquifer.bottom[j])
    end
    return saturation * aquifer.conductance[nzi]
end

function flux!(Q, aquifer, connectivity)
    for i in 1:connectivity.ncell
        # Loop over connections for cell j
        for nzi in connections(connectivity, i)
            # connection from i -> j
            j = connectivity.rowval[nzi]
            Δϕ = aquifer.head[i] - aquifer.head[j]
            cond = conductance(aquifer, i, j, nzi)
            Q[i] -= cond * Δϕ
        end
    end
    return Q
end

struct ConstantHead{T}
    head::Vector{T}
    index::Vector{Int}
end

"""
    stable_timestep(aquifer)

Compute a stable timestep size given the forward-in-time, central in space scheme.
The following criterion can be found in Chu & Willis (1984)

Δt * k * H / (Δx * Δy * S) <= 1/4
"""
function stable_timestep(aquifer)
    Δtₘᵢₙ = Inf
    for i in eachindex(aquifer.head)
        Δt = aquifer.area[i] * storativity(aquifer)[i] / (aquifer.k[i] * saturated_thickness(aquifer, i))
        Δtₘᵢₙ = Δt < Δtₘᵢₙ ? Δt : Δtₘᵢₙ
    end
    return 0.25 * Δtₘᵢₙ
end

minimum_head(aquifer::ConfinedAquifer) = aquifer.head
minimum_head(aquifer::UnconfinedAquifer) = max.(aquifer.head, aquifer.bottom)

function update(gwf, Q, Δt)
    Q .= 0.0  # TODO: Probably remove this when linking with other components
    flux!(Q, gwf.aquifer, gwf.connectivity)
    for boundary in gwf.boundaries
        flux!(Q, boundary, gwf.aquifer)
    end
    gwf.aquifer.head .+= (Q ./ gwf.aquifer.area .* Δt ./ storativity(gwf.aquifer))
    # Set constant head (dirichlet) boundaries
    gwf.aquifer.head[gwf.constanthead.index] .= gwf.constanthead.head
    # Make sure no heads ends up below an unconfined aquifer bottom
    gwf.aquifer.head .= minimum_head(gwf.aquifer)
    return gwf
end

Base.@kwdef struct GroundwaterFlow{A, B}
    aquifer::A
    connectivity::Connectivity
    constanthead::ConstantHead
    boundaries::Vector{B}
    function GroundwaterFlow(
            aquifer::A,
            connectivity,
            constanthead,
            boundaries::Vector{B}) where {A <: Aquifer, B <: AquiferBoundaryCondition}
        initialize_conductance!(aquifer, connectivity)
        new{A, B}(aquifer, connectivity, constanthead, boundaries)
    end
end