text:
#=
-------------------------------------------------------------------------------------------------------------
This file contains the student template for Project 6: 2D Wave Equation using the unified CG/DG method with AGGP storage
introduced in F.X. Giraldo's Introduction to Element-based Galerkin Methods using
Tensor-Product Bases: Analysis, Algorithms, and Applications.

The solution is obtained using Algorithm 17.7 whereby the RHS vector is constructed at each time stage.
The volume integral contribution is constructed using Algorithm 16.11 except that the element differentiation matrix is computed and stored.
The flux integral contribution is constructed using Algorithm 16.12 (all in weak form).

It uses the LSRK45 method as the time-integrator and Rusanov fluxes for DG.

Written by F.X. Giraldo on July 10, 2021.
Department of Applied Mathematics
Naval Postgraduate School
Monterey, CA 93943

The interpolation points used are the following:
ipoints=1: Lobatto
ipoints=2: Legendre
ipoints=3: Chebyshev
ipoints=4: Equi-spaced

The integration points used are the following:
qpoints=1: Lobatto
qpoints=2: Legendre
-------------------------------------------------------------------------------------------------------------
=#

using Plots
using LinearAlgebra

include("QuadraturePoints.jl")
include("element_matrices.jl")
include("global_matrices.jl")
include("create_grid.jl")
include("create_CGDG_storage.jl")
include("create_side.jl")
include("create_face.jl")
include("compute_normals.jl")
include("create_face_periodicity.jl")
include("compute_metrics.jl")
include("global_matrices.jl")
include("exact_solution.jl")
include("compute_norm.jl")
include("ti_LSRK.jl")

#Some Constants
DFloat = Float64
machine_zero=eps(DFloat)

#Main Function
function main()

    #------------------------Only Change stuff here-------------------#
    N=8 #polynomial order
    Nex=4 #number of elements along x (ξ)
    Ney=4 #number of elements along y (η)
    space_method="DG" #CG or DG
    time_final=DFloat(0.25) #in revolutions
    plot_grid=false
    plot_solution=true
    plot_movie=false
    warp_grid=false
    #------------------------Only Change stuff here-------------------#

    Q=N #changing this requires modifying CREATE_RHS_FLUX which assumes Q=N
    Np=N+1
    Nq=Q+1
    Npts=Np^2
    Nqs=Nq^2
    ipoints=1 #lobatto
    qpoints=1 #lobatto
    Courant_max=DFloat(0.5)
    icase=1 #CCW moving Gaussian
    Ne=Nex*Ney
    Nx=Nex*N+1
    Ny=Ney*N+1
    Npoin_CG=Nx*Ny
    Npoin_DG=Npts*Ne
    Npoin=Npoin_CG
    Nboun=2*Nex + 2*Ney
    Nface=2*Ne + Nex + Ney
    if icase == 1
        c=2*π;
    end
    rotate_grid=false #not working with periodicity
    time_final=time_final*c
    @show(N,Q,Ne,Npoin_CG,Npoin_DG,Nboun,Nface,Nx,Ny)

    #Select Interpolation Points
    if ipoints == 1
        (ξ,ω) = QuadraturePoints.lobatto_gauss(Np,DFloat) #Lobatto
    elseif ipoints == 2
        (ξ,ω) = QuadraturePoints.legendre_gauss(Np,DFloat) #Legendre
    elseif ipoints == 3
        (ξ,ω) = QuadraturePoints.chebyshev_gauss(Np,DFloat) #Chebyshev
    elseif ipoints == 4
        (ξ,ω) = QuadraturePoints.equispaced_points(Np,DFloat) #Equi-spaced
    end

    #Select Quadrature Points
    if qpoints == 1
        (ξq,ωq) = QuadraturePoints.lobatto_gauss(Nq,DFloat) #Lobatto
    elseif qpoints == 2
        (ξq,ωq) = QuadraturePoints.legendre_gauss(Nq,DFloat) #Legendre
    end

    #Construct Lagrange Polynomials
    (ψ,dψ) = QuadraturePoints.lagrange_basis(Np,Nq,ξ,ξq,DFloat)

    #Construct Grid and (i,e)->I map
    (coord_CG,intma_CG,bsido_CG,periodicity_CG) = create_grid(Np,Npoin_CG,Ne,Nboun,Nex,Ney,ξ,plot_grid,warp_grid,rotate_grid,DFloat)

    #Store the CGDG-Storage Grid
    (coord,intma,periodicity,DG_to_CG,Npoin) = create_CGDG_Storage(space_method,Npoin_CG,Npoin_DG,coord_CG,intma_CG,periodicity_CG,Ne,Np,DFloat)
    @show(space_method,Npoin)

    #Construct Metric Terms
    (ξ_x,ξ_y,η_x,η_y,jac) = compute_metrics(coord_CG,intma_CG,ψ,dψ,Ne,Np,Nq,DFloat)

    #Construct Face arrays
    (iside,~) = create_side(intma_CG,bsido_CG,Npoin_CG,Ne,Nboun,Nface,Np,DFloat)
    (face,mapL,mapR) = create_face(iside,intma_CG,Nface,Np)
    (normals,jac_face) = compute_normals(face,intma_CG,coord_CG,Nface,Np,Nq,ψ,dψ)
    (face,iside) = create_face_periodicity!(face,iside,coord_CG,Nface,Nboun)

    #Construct Element Matrices
    (Me,De_x,De_y) = element_matrices(ψ,dψ,ξ_x,ξ_y,η_x,η_y,jac,ωq,Ne,Np,Nq,DFloat)

    #Construct Global Matrices
    M = global_matrices(Me,intma,periodicity,Ne,Np,Npoin,DFloat)

    #Compute Exact Solution
    time=DFloat(0.0)
    (qe,ue) = exact_solution(coord,Npoin,time,icase,DFloat)
    q0=copy(qe)

    #Recompute DT, NTIME
    dx_min=coord[1,2]-coord[1,1]
    umax=maximum(ue)
    dt=Courant_max*dx_min/umax
    courant=umax*dt/dx_min
    ntime = ceil(Int64, time_final / dt)
    println("Courant = ",courant," dt = ",dt," ntime = ",ntime)

    #Time Integration
    #-----------------------------Students Add their Functions inside ti_LSRK--------------------------#
    #-----------------------------Students Add their Functions inside ti_LSRK--------------------------#
    (q0,time) = ti_LSRK!(q0,ue,coord,M,De_x,De_y,intma,periodicity,face,mapL,mapR,normals,jac_face,ωq,time,ntime,dt,space_method,plot_movie,Nx,Ny,coord_CG,periodicity_CG,DG_to_CG,DFloat)
    #-----------------------------Students Add their Functions inside ti_LSRK--------------------------#
    #-----------------------------Students Add their Functions inside ti_LSRK--------------------------#

    #Print Time and Extrema
    println(" itime = ",ntime," time = ",time," qmax = ",maximum(q0)," qmin = ",minimum(q0))

    #Compute L2 Norm
    (qe,ue) = exact_solution(coord,Npoin,time,icase,DFloat)
    norm2 = compute_norm(q0,qe,Npoin,DFloat)
    println("norm.L2 = ",norm2)

    #Plot Solution
    if (plot_solution)
        #Compute gridpoint solution
        q_sol=zeros(DFloat,Npoin_CG)
        qe_sol=zeros(DFloat,Npoin_CG)
        lhowm=zeros(Int64,Npoin_CG)
        for i=1:Npoin
            ip_CG=periodicity_CG[DG_to_CG[i]]
            q_sol[ip_CG] += q0[i]
            qe_sol[ip_CG] += qe[i]
            lhowm[ip_CG] += 1
        end
        for i=1:Npoin_CG
            q_sol[i]=q_sol[i]/lhowm[i]
            qe_sol[i]=qe_sol[i]/lhowm[i];
        end

        x=reshape(coord_CG[1,:],(Nx,Ny))
        y=reshape(coord_CG[2,:],(Nx,Ny))
        xx=x[:,1]
        yy=y[1,:]

        #Numerical Solution
        qq=reshape(q_sol,(Nx,Ny))
        p1 = contourf(xx,yy,qq')
        title!(p1,"$space_method solution")
        xlabel!(p1,"x")
        ylabel!(p1,"y")

        #Exact Solution
        qq=reshape(qe_sol,(Nx,Ny))
        p2 = contourf(xx,yy,qq')
        title!(p2,"Exact solution")
        xlabel!(p2,"x")
        ylabel!(p2,"y")
        figure_1=plot(p1,p2)
        display(figure_1)
    end

    println("Done") #output

end

#----------------------------------#
# Run the main function
#----------------------------------#
main()


text:
"""
    lagrange_point(sp::AbstractVector{T}, x) where {T<:Real}

Calculate Legendre polynomials of solution points sp at location x
"""
function lagrange_point(sp::AbstractVector{T}, x::Real) where {T<:Real}
    l = similar(sp)
    nsp = length(sp)

    for k = 1:nsp
        tmp = 1.0
        for j = 1:nsp
            if j != k
                tmp *= (x - sp[j]) / (sp[k] - sp[j])
            end
        end
        l[k] = tmp
    end

    return l
end

function lagrange_point(sp, x::AbstractVector{T}) where {T<:Real}
    lp = zeros(eltype(sp), axes(x, 1), axes(sp, 1))

    for i in axes(lp, 1)
        lp[i, :] .= lagrange_point(sp, x[i])
    end

    return lp
end

"""
    ∂lagrange(sp::T) where {T<:AbstractVector{<:Real}}
    ∂lagrange(V, Vr)
    ∂lagrange(V, Vr, Vs)

Calculate derivatives of Lagrange polynomials dlⱼ(rᵢ)
"""
function ∂lagrange(sp::T) where {T<:AbstractVector{<:Real}}
    nsp = length(sp)
    lpdm = similar(sp, nsp, nsp)

    for k = 1:nsp, m = 1:nsp
        lsum = 0.0
        for l = 1:nsp
            tmp = 1.0
            for j = 1:nsp
                if j != k && j != l
                    tmp *= (sp[m] - sp[j]) / (sp[k] - sp[j])
                end
            end
            if l != k
                lsum += tmp / (sp[k] - sp[l])
            end
        end
        lpdm[m, k] = lsum
    end

    return lpdm
end

# ------------------------------------------------------------
# Vandermonde matrix based evaluation
# ------------------------------------------------------------

function ∂lagrange(V, Vr)
    Np = size(V, 1)

    ∂l = zeros(Np, Np)
    for i = 1:Np
        ∂l[i, :] .= V' \ Vr[i, :]
    end

    return ∂l
end

function ∂lagrange(V, Vr, Vs)
    Np = size(V, 1)

    ∂l = zeros(Np, Np, 2)
    for i = 1:Np
        ∂l[i, :, 1] .= V' \ Vr[i, :]
        ∂l[i, :, 2] .= V' \ Vs[i, :]
    end

    return ∂l
end

"""
One-shot calculation of derivatives of Lagrange polynomials and the values at cell edge

"""
function standard_lagrange(x)
    ll = lagrange_point(x, -1.0)
    lr = lagrange_point(x, 1.0)
    lpdm = ∂lagrange(x)

    return ll, lr, lpdm
end

ext: jl


text:
###############################################################################
#
#   fmpq_poly.jl : Flint polynomials over fmpq
#
###############################################################################

export FmpqPolyRing, fmpq_poly

###############################################################################
#
#   Data type and parent methods
#
###############################################################################

parent_type(::Type{fmpq_poly}) = FmpqPolyRing

elem_type(::Type{FmpqPolyRing}) = fmpq_poly

base_ring(a::FmpqPolyRing) = a.base_ring

var(a::FmpqPolyRing) = a.S

###############################################################################
#
#   Basic manipulation
#
###############################################################################

doc"""
    denominator(a::fmpq_poly)
> Return the least common denominator of the coefficients of the polynomial
> $a$.
"""
function denominator(a::fmpq_poly)
   z = fmpz()
   ccall((:fmpq_poly_get_denominator, :libflint), Void,
         (Ref{fmpz}, Ref{fmpq_poly}), z, a)
   return z
end

length(x::fmpq_poly) = ccall((:fmpq_poly_length, :libflint), Int,
                                   (Ref{fmpq_poly},), x)

set_length!(x::fmpq_poly, n::Int) = ccall((:_fmpq_poly_set_length, :libflint), Void,
                                   (Ref{fmpq_poly}, Int), x, n)

function coeff(x::fmpq_poly, n::Int)
   n < 0 && throw(DomainError())
   z = fmpq()
   ccall((:fmpq_poly_get_coeff_fmpq, :libflint), Void,
               (Ref{fmpq}, Ref{fmpq_poly}, Int), z, x, n)
   return z
end

zero(a::FmpqPolyRing) = a(0)

one(a::FmpqPolyRing) = a(1)

gen(a::FmpqPolyRing) = a([zero(base_ring(a)), one(base_ring(a))])

isgen(x::fmpq_poly) = ccall((:fmpq_poly_is_x, :libflint), Bool,
                            (Ref{fmpq_poly},), x)

function deepcopy_internal(a::fmpq_poly, dict::ObjectIdDict)
   z = fmpq_poly(a)
   z.parent = parent(a)
   return z
end

###############################################################################
#
#   Canonicalisation
#
###############################################################################

canonical_unit(a::fmpq_poly) = canonical_unit(lead(a))

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

#=
#use the generic one to print //
function show(io::IO, x::fmpq_poly)
   if length(x) == 0
      print(io, "0")
   else
      cstr = ccall((:fmpq_poly_get_str_pretty, :libflint), Ref{UInt8},
          (Ref{fmpq_poly}, Ref{UInt8}), x, string(var(parent(x))))

      print(io, unsafe_string(cstr))

      ccall((:flint_free, :libflint), Void, (Ref{UInt8},), cstr)
   end
end
=#

function show(io::IO, p::FmpqPolyRing)
   print(io, "Univariate Polynomial Ring in ")
   print(io, string(var(p)))
   print(io, " over ")
   show(io, p.base_ring)
end

show_minus_one(::Type{fmpq_poly}) = show_minus_one(FracElem{fmpz})

###############################################################################
#
#   Unary operations
#
###############################################################################

function -(x::fmpq_poly)
   z = parent(x)()
   ccall((:fmpq_poly_neg, :libflint), Void,
         (Ref{fmpq_poly}, Ref{fmpq_poly}), z, x)
   return z
end

###############################################################################
#
#   Binary operations
#
###############################################################################

function +(x::fmpq_poly, y::fmpq_poly)
   check_parent(x, y)
   z = parent(x)()
   ccall((:fmpq_poly_add, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly},  Ref{fmpq_poly}),
               z, x, y)
   return z
end

function -(x::fmpq_poly, y::fmpq_poly)
   check_parent(x, y)
   z = parent(x)()
   ccall((:fmpq_poly_sub, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly},  Ref{fmpq_poly}),
               z, x, y)
   return z
end

function *(x::fmpq_poly, y::fmpq_poly)
   check_parent(x, y)
   z = parent(x)()
   ccall((:fmpq_poly_mul, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly},  Ref{fmpq_poly}),
               z, x, y)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::Int, y::fmpq_poly)
   z = parent(y)()
   ccall((:fmpq_poly_scalar_mul_si, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Int), z, y, x)
   return z
end

function *(x::fmpz, y::fmpq_poly)
   z = parent(y)()
   ccall((:fmpq_poly_scalar_mul_fmpz, :libflint), Void,
         (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpz}), z, y, x)
   return z
end

function *(x::fmpq, y::fmpq_poly)
   z = parent(y)()
   ccall((:fmpq_poly_scalar_mul_fmpq, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq}), z, y, x)
   return z
end

function +(x::fmpq_poly, y::Int)
   z = parent(x)()
   ccall((:fmpq_poly_add_si, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Int), z, x, y)
   return z
end

function +(x::fmpq_poly, y::fmpz)
   z = parent(x)()
   ccall((:fmpq_poly_add_fmpz, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpz}), z, x, y)
   return z
end

function +(x::fmpq_poly, y::fmpq)
   z = parent(x)()
   ccall((:fmpq_poly_add_fmpq, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq}), z, x, y)
   return z
end

function -(x::fmpq_poly, y::Int)
   z = parent(x)()
   ccall((:fmpq_poly_sub_si, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Int), z, x, y)
   return z
end

function -(x::fmpq_poly, y::fmpz)
   z = parent(x)()
   ccall((:fmpq_poly_sub_fmpz, :libflint), Void,
         (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpz}), z, x, y)
   return z
end

function -(x::fmpq_poly, y::fmpq)
   z = parent(x)()
   ccall((:fmpq_poly_sub_fmpq, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq}), z, x, y)
   return z
end

function -(x::Int, y::fmpq_poly)
   z = parent(y)()
   ccall((:fmpq_poly_si_sub, :libflint), Void,
                (Ref{fmpq_poly}, Int, Ref{fmpq_poly}), z, x, y)
   return z
end

function -(x::fmpz, y::fmpq_poly)
   z = parent(y)()
   ccall((:fmpq_poly_fmpz_sub, :libflint), Void,
         (Ref{fmpq_poly}, Ref{fmpz}, Ref{fmpq_poly}), z, x, y)
   return z
end

function -(x::fmpq, y::fmpq_poly)
   z = parent(y)()
   ccall((:fmpq_poly_fmpq_sub, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq}, Ref{fmpq_poly}), z, x, y)
   return z
end

+(x::Int, y::fmpq_poly) = y + x

+(x::fmpz, y::fmpq_poly) = y + x

+(x::fmpq, y::fmpq_poly) = y + x

*(x::fmpq_poly, y::Int) = y*x

*(x::fmpq_poly, y::fmpz) = y*x

*(x::fmpq_poly, y::fmpq) = y*x

+(x::Integer, y::fmpq_poly) = y + fmpz(x)

-(x::Integer, y::fmpq_poly) = fmpz(x) - y

*(x::Integer, y::fmpq_poly) = fmpz(x)*y

+(x::fmpq_poly, y::Integer) = x + fmpz(y)

-(x::fmpq_poly, y::Integer) = x - fmpz(y)

*(x::fmpq_poly, y::Integer) = fmpz(y)*x

+(x::Rational, y::fmpq_poly) = fmpq(x) + y

-(x::Rational, y::fmpq_poly) = fmpq(x) - y

*(x::Rational, y::fmpq_poly) = fmpq(x) * y

+(x::fmpq_poly, y::Rational) = x + fmpq(y)

-(x::fmpq_poly, y::Rational) = x - fmpq(y)

*(x::fmpq_poly, y::Rational) = x * fmpq(y)

###############################################################################
#
#   Powering
#
###############################################################################

function ^(x::fmpq_poly, y::Int)
   y < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fmpq_poly_pow, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Int),
               z, x, y)
   return z
end

###############################################################################
#
#   Comparisons
#
###############################################################################

function ==(x::fmpq_poly, y::fmpq_poly)
   check_parent(x, y)
   return ccall((:fmpq_poly_equal, :libflint), Bool,
                                      (Ref{fmpq_poly}, Ref{fmpq_poly}), x, y)
end

###############################################################################
#
#   Ad hoc comparisons
#
###############################################################################

function ==(x::fmpq_poly, y::fmpq)
   if length(x) > 1
      return false
   elseif length(x) == 1
      z = fmpq()
      ccall((:fmpq_poly_get_coeff_fmpq, :libflint), Void,
                       (Ref{fmpq}, Ref{fmpq_poly}, Int), z, x, 0)
      return ccall((:fmpq_equal, :libflint), Bool,
               (Ref{fmpq}, Ref{fmpq}, Int), z, y, 0)
   else
      return iszero(y)
   end
end

==(x::fmpq, y::fmpq_poly) = y == x

==(x::fmpq_poly, y::Rational{T}) where T <: Union{Int, BigInt} = x == fmpq(y)

==(x::Rational{T}, y::fmpq_poly) where T <: Union{Int, BigInt} = y == x

###############################################################################
#
#   Truncation
#
###############################################################################

function truncate(a::fmpq_poly, n::Int)
   n < 0 && throw(DomainError())

   if length(a) <= n
      return a
   end

   z = parent(a)()
   ccall((:fmpq_poly_set_trunc, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Int), z, a, n)
   return z
end

function mullow(x::fmpq_poly, y::fmpq_poly, n::Int)
   check_parent(x, y)
   n < 0 && throw(DomainError())

   z = parent(x)()
   ccall((:fmpq_poly_mullow, :libflint), Void,
         (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly}, Int), z, x, y, n)
   return z
end

###############################################################################
#
#   Reversal
#
###############################################################################

function reverse(x::fmpq_poly, len::Int)
   len < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fmpq_poly_reverse, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Int), z, x, len)
   return z
end

###############################################################################
#
#   Shifting
#
###############################################################################

function shift_left(x::fmpq_poly, len::Int)
   len < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fmpq_poly_shift_left, :libflint), Void,
      (Ref{fmpq_poly}, Ref{fmpq_poly}, Int), z, x, len)
   return z
end

function shift_right(x::fmpq_poly, len::Int)
   len < 0 && throw(DomainError())
   z = parent(x)()
   ccall((:fmpq_poly_shift_right, :libflint), Void,
       (Ref{fmpq_poly}, Ref{fmpq_poly}, Int), z, x, len)
   return z
end

###############################################################################
#
#   Euclidean division
#
###############################################################################

function mod(x::fmpq_poly, y::fmpq_poly)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   r = parent(x)()
   ccall((:fmpq_poly_rem, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly}),
               r, x, y)
   return r
end

rem(x::fmpq_poly, y::fmpq_poly) = mod(x, y)

function divrem(x::fmpq_poly, y::fmpq_poly)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   q = parent(x)()
   r = parent(x)()
   ccall((:fmpq_poly_divrem, :libflint), Void,
         (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly}),
               q, r, x, y)
   return q, r
end

###############################################################################
#
#   Exact division
#
###############################################################################

function div(x::fmpq_poly, y::fmpq_poly)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   ccall((:fmpq_poly_div, :libflint), Void,
            (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly}), z, x, y)
   return z
end

divexact(x::fmpq_poly, y::fmpq_poly) = div(x,y)

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::fmpq_poly, y::fmpz)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   ccall((:fmpq_poly_scalar_div_fmpz, :libflint), Void,
          (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpz}), z, x, y)
   return z
end

function divexact(x::fmpq_poly, y::fmpq)
   iszero(y) && throw(DivideError())
   z = parent(x)()
   ccall((:fmpq_poly_scalar_div_fmpq, :libflint), Void,
          (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq}), z, x, y)
   return z
end

function divexact(x::fmpq_poly, y::Int)
   y == 0 && throw(DivideError())
   z = parent(x)()
   ccall((:fmpq_poly_scalar_div_si, :libflint), Void,
                        (Ref{fmpq_poly}, Ref{fmpq_poly}, Int), z, x, y)
   return z
end

divexact(x::fmpq_poly, y::Integer) = divexact(x, fmpz(y))

divexact(x::fmpq_poly, y::Rational{T}) where T <: Union{Int, BigInt} = divexact(x, fmpq(y))

###############################################################################
#
#   Removal and valuation
#
###############################################################################

function divides(z::fmpq_poly, x::fmpq_poly)
   q, r = divrem(z, x)
   return iszero(r), q
end

###############################################################################
#
#   Content, primitive part, GCD and LCM
#
###############################################################################

function gcd(x::fmpq_poly, y::fmpq_poly)
   check_parent(x, y)
   z = parent(x)()
   ccall((:fmpq_poly_gcd, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly}), z, x, y)
   return z
end

function content(x::fmpq_poly)
   z = fmpq()
   ccall((:fmpq_poly_content, :libflint), Void,
         (Ref{fmpq}, Ref{fmpq_poly}), z, x)
   return z
end

function primpart(x::fmpq_poly)
   z = parent(x)()
   ccall((:fmpq_poly_primitive_part, :libflint), Void,
         (Ref{fmpq_poly}, Ref{fmpq_poly}), z, x)
   return z
end

###############################################################################
#
#   Evaluation
#
###############################################################################

function evaluate(x::fmpq_poly, y::fmpz)
   z = fmpq()
   ccall((:fmpq_poly_evaluate_fmpz, :libflint), Void,
                (Ref{fmpq}, Ref{fmpq_poly}, Ref{fmpz}), z, x, y)
   return z
end

function evaluate(x::fmpq_poly, y::fmpq)
   z = fmpq()
   ccall((:fmpq_poly_evaluate_fmpq, :libflint), Void,
                (Ref{fmpq}, Ref{fmpq_poly}, Ref{fmpq}), z, x, y)
   return z
end

evaluate(x::fmpq_poly, y::Integer) = evaluate(x, fmpz(y))

evaluate(x::fmpq_poly, y::Rational) = evaluate(x, fmpq(y))

###############################################################################
#
#   Composition
#
###############################################################################

function compose(x::fmpq_poly, y::fmpq_poly)
   check_parent(x, y)
   z = parent(x)()
   ccall((:fmpq_poly_compose, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly}), z, x, y)
   return z
end

###############################################################################
#
#   Derivative
#
###############################################################################

function derivative(x::fmpq_poly)
   z = parent(x)()
   ccall((:fmpq_poly_derivative, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}), z, x)
   return z
end

###############################################################################
#
#   Integral
#
###############################################################################

function integral(x::fmpq_poly)
   z = parent(x)()
   ccall((:fmpq_poly_integral, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}), z, x)
   return z
end

###############################################################################
#
#   Resultant
#
###############################################################################

function resultant(x::fmpq_poly, y::fmpq_poly)
   check_parent(x, y)
   z = fmpq()
   ccall((:fmpq_poly_resultant, :libflint), Void,
                (Ref{fmpq}, Ref{fmpq_poly}, Ref{fmpq_poly}), z, x, y)
   return z
end

###############################################################################
#
#   GCDX
#
###############################################################################

function gcdx(x::fmpq_poly, y::fmpq_poly)
   check_parent(x, y)
   z = parent(x)()
   u = parent(x)()
   v = parent(x)()
   ccall((:fmpq_poly_xgcd, :libflint), Void,
        (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly},
                                     Ref{fmpq_poly}), z, u, v, x, y)
   return (z, u, v)
end

################################################################################
#
#   Factorization
#
################################################################################

doc"""
    factor(x::fmpq_poly)
> Returns the factorization of $x$.
"""
function factor(x::fmpq_poly)
   res, z = _factor(x)
   return Fac(parent(x)(z), res)
end

function _factor(x::fmpq_poly)
   res = Dict{fmpq_poly, Int}()
   y = fmpz_poly()
   ccall((:fmpq_poly_get_numerator, :libflint), Void,
         (Ref{fmpz_poly}, Ref{fmpq_poly}), y, x)
   fac = fmpz_poly_factor()
   ccall((:fmpz_poly_factor, :libflint), Void,
              (Ref{fmpz_poly_factor}, Ref{fmpz_poly}), fac, y)
   z = fmpz()
   ccall((:fmpz_poly_factor_get_fmpz, :libflint), Void,
            (Ref{fmpz}, Ref{fmpz_poly_factor}), z, fac)
   f = fmpz_poly()
   for i in 1:fac.num
      ccall((:fmpz_poly_factor_get_fmpz_poly, :libflint), Void,
            (Ref{fmpz_poly}, Ref{fmpz_poly_factor}, Int), f, fac, i - 1)
      e = unsafe_load(fac.exp, i)
      res[parent(x)(f)] = e
   end
   return res, fmpq(z, denominator(x))
end

###############################################################################
#
#   Signature
#
###############################################################################

doc"""
    signature(f::fmpq_poly)
> Return the signature of $f$, i.e. a tuple $(r, s)$ where $r$ is the number of
> real roots of $f$ and $s$ is half the number of complex roots.
"""
function signature(f::fmpq_poly)
   r = Array{Int}(1)
   s = Array{Int}(1)
   z = fmpz_poly()
   ccall((:fmpq_poly_get_numerator, :libflint), Void,
         (Ref{fmpz_poly}, Ref{fmpq_poly}), z, f)
   return signature(z)
end

###############################################################################
#
#   Speedups for polynomials over fmpq_polys
#
###############################################################################

function *(a::Generic.Poly{fmpq_poly}, b::Generic.Poly{fmpq_poly})
   check_parent(a, b)
   if min(length(a), length(b)) < 40
      return mul_classical(a, b)
   else
      return mul_ks(a, b)
   end
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(z::fmpq_poly)
   ccall((:fmpq_poly_zero, :libflint), Void,
                    (Ref{fmpq_poly},), z)
   return z
end

function fit!(z::fmpq_poly, n::Int)
   ccall((:fmpq_poly_fit_length, :libflint), Void,
                    (Ref{fmpq_poly}, Int), z, n)
   return nothing
end

function setcoeff!(z::fmpq_poly, n::Int, x::fmpz)
   ccall((:fmpq_poly_set_coeff_fmpz, :libflint), Void,
                    (Ref{fmpq_poly}, Int, Ref{fmpz}), z, n, x)
   return z
end

function setcoeff!(z::fmpq_poly, n::Int, x::fmpq)
   ccall((:fmpq_poly_set_coeff_fmpq, :libflint), Void,
                    (Ref{fmpq_poly}, Int, Ref{fmpq}), z, n, x)
   return z
end

function mul!(z::fmpq_poly, x::fmpq_poly, y::fmpq_poly)
   ccall((:fmpq_poly_mul, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly}), z, x, y)
   return z
end

function addeq!(z::fmpq_poly, x::fmpq_poly)
   ccall((:fmpq_poly_add, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly}), z, z, x)
   return z
end

function add!(z::fmpq_poly, x::fmpq_poly, y::fmpq_poly)
   ccall((:fmpq_poly_add, :libflint), Void,
                (Ref{fmpq_poly}, Ref{fmpq_poly}, Ref{fmpq_poly}), z, x, y)
   return z
end

###############################################################################
#
#   Promotions
#
###############################################################################

promote_rule(::Type{fmpq_poly}, ::Type{T}) where {T <: Integer} = fmpq_poly

promote_rule(::Type{fmpq_poly}, ::Type{fmpz}) = fmpq_poly

promote_rule(::Type{fmpq_poly}, ::Type{fmpq}) = fmpq_poly

promote_rule(::Type{fmpq_poly}, ::Type{Rational{T}}) where T <: Union{Int, BigInt} = fmpq_poly

###############################################################################
#
#   Polynomial substitution
#
###############################################################################

(f::fmpq_poly)(a::fmpq) = evaluate(f, a)

(f::fmpq_poly)(a::Rational) = evaluate(f, fmpq(a))

###############################################################################
#
#   Parent object call overloads
#
###############################################################################

function (a::FmpqPolyRing)()
   z = fmpq_poly()
   z.parent = a
   return z
end

function (a::FmpqPolyRing)(b::Int)
   z = fmpq_poly(b)
   z.parent = a
   return z
end

function (a::FmpqPolyRing)(b::Integer)
   z = fmpq_poly(fmpz(b))
   z.parent = a
   return z
end

function (a::FmpqPolyRing)(b::fmpz)
   z = fmpq_poly(b)
   z.parent = a
   return z
end

function (a::FmpqPolyRing)(b::fmpq)
   z = fmpq_poly(b)
   z.parent = a
   return z
end

function (a::FmpqPolyRing)(b::Array{fmpq, 1})
   z = fmpq_poly(b)
   z.parent = a
   return z
end

(a::FmpqPolyRing)(b::Rational) = a(fmpq(b))

(a::FmpqPolyRing)(b::Array{T, 1}, copy::Bool=true) where {T <: Integer} = a(map(fmpq, b))

(a::FmpqPolyRing)(b::Array{Rational{T}, 1}, copy::Bool=true) where {T <: Integer} = a(map(fmpq, b))

(a::FmpqPolyRing)(b::Array{fmpz, 1}, copy::Bool=true) = a(map(fmpq, b))

(a::FmpqPolyRing)(b::fmpq_poly) = b

function (a::FmpqPolyRing)(b::fmpz_poly)
   z = fmpq_poly(b)
   z.parent = a
   return z
end

###############################################################################
#
#   PolynomialRing constructor
#
###############################################################################

function PolynomialRing(R::FlintRationalField, s::AbstractString; cached = true)
   S = Symbol(s)

   parent_obj = FmpqPolyRing(R, S, cached)

   return parent_obj, parent_obj([fmpq(0), fmpq(1)])
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE

# Kaltenbacher, Manfred. Numerical simulation of mechatronic sensors and actuators: finite elements for computational multiphysics. Springer, 2015.
code = create_basis(
    :Pyr5A,
    "5 node linear pyramid element",
    (
     ( 1.0,  1.0, 0.0), # N1
     ( 1.0, -1.0, 0.0), # N2
     (-1.0, -1.0, 0.0), # N3
     (-1.0,  1.0, 0.0), # N4
     ( 0.0,  0.0, 1.0), # N5
    ),
    (
     "1/4*( (1+u)*(1+v) - w + u*v*w/(1-w) )",
     "1/4*( (1+u)*(1-v) - w + u*v*w/(1-w) )",
     "1/4*( (1-u)*(1-v) - w + u*v*w/(1-w) )",
     "1/4*( (1-u)*(1+v) - w + u*v*w/(1-w) )",
     "w",
    ),
   )
eval(code)

# source: Code Aster documentation?
code = create_basis(
    :Pyr5,
    "5 node linear pyramid element",
    (
     (-1.0, -1.0, -1.0), # N1
     ( 1.0, -1.0, -1.0), # N2
     ( 1.0,  1.0, -1.0), # N3
     (-1.0,  1.0, -1.0), # N4
     ( 0.0,  0.0,  1.0), # N5
    ),
    (
     "1/8 * (1-u) * (1-v) * (1-w)",
     "1/8 * (1+u) * (1-v) * (1-w)",
     "1/8 * (1+u) * (1+v) * (1-w)",
     "1/8 * (1-u) * (1+v) * (1-w)",
     "1/2 * (1+w)",
    ),
   )
eval(code)


text:
function writevtk(X::LazySet{N}; file="output") where {N}
    points, connec = triangulate(X)
    ntriangles = length(connec)
    cell = VTKPolyhedron(1:ntriangles, connec...)

    vtk_grid(file, points, [cell]; compress=false) do vtk
        #
    end
end

function writevtk(X::AbstractVector{VT}; file="output") where {N, VT<:LazySet{N}}
    points_list = Vector{Matrix{Float32}}()
    cell_list = Vector{VTKPolyhedron}()
    count = 0
    for Xi in X
        points, connec = triangulate(Xi)
        ntriangles = length(connec)
        connec = [c .+ count for c in connec]
        cell = VTKPolyhedron(1:ntriangles, connec...)
        count += size(points, 2)
        push!(points_list, points)
        push!(cell_list, cell)
    end
    points_matrix = reduce(hcat, points_list)

    vtk_grid(file, points_matrix,  cell_list; compress=false) do vtk
        #
    end
end


text:
using Plots

function (D(R))
    D = (3*G*M*mdot)/(8*pi*R^3)*(1-(R_isco/R)^(1/2))
    return D
end

G = 6.67e-11
M = 2e30*10
c = 3e8
eta = 0.1

R_isco = 6*G*M/(c^2)

L_edd = 3e4 *3.8e26 * (2e30*10)/1.99e30

Mdot = L_edd/(c^2*eta)
mdot = 0.1*Mdot

r_vals = LinRange(R_isco,10*R_isco,10000)
D_vals = []

for R in r_vals #creates a series of D (dissipation) values corresponding to radial distances
    push!(D_vals, D(R)) #(cycles through R values, taken from r_vals,
                        # and inputs them into function D(R), which is a function of R
                        # the resulting value is then added to D_vals)
end

sig_SB = 5.67e-8

T_vals = []
function (T(D))
    T = (D/sig_SB)^(1/4)
    return T
end

for D in D_vals
    push!(T_vals, T(D))

end

plot(r_vals, T_vals)


text:
# This includes routines that deal with tracers

"""
    Structure that has information about the tracers

    General form:
        Tracer(num, coord, T)

        with:

            num:   number of the tracer (integer)

            coord: coordinates of the tracer
                    2D - [x; z]
                    3D - [x; y; z]

            T:      Temperature of the tracer [Celcius]
"""
@with_kw struct Tracer
    num         ::  Int64     =  0           # number
    coord       ::  Vector{Float64}          # holds coordinates [2D or 3D]
    T           ::  Float64   =  900         # temperature
    Phase       ::  Int64     =  1           # Phase (aka rock type) of the Tracer
    Phi_melt    ::  Float64   =  0           # Melt fraction on Tracers
end
#    Chemistry   ::  Vector{Float64} = []    # Could @ some stage hold the evolving chemistry of the magma

"""

    Function that updates properties on tracers

        General form:

        Tracers = UpdateTracers(Tracers, Grid, T, Phi, InterpolationMethod);

        with:
            Tracers:   StructArray that contains tracers

            Grid:   regular grid on which the parameters to be interpolated are defined
                    2D - (X,Z)
                    3D - (X,Y,Z)

            T:      Temperature that is defined on the grid.

            Phi:    Solid fraction defined on grid

            InterpolationMethod:    Interpolation method from grid->Tracers
                    "Cubic"     -   Cubic interpolation
                    "Quadratic" -   Quadratic interpolation (default)
                    "Linear"    -   Linear interpolation

        out:
            Tracers:    Tracers structure with updated T and melt fraction fields
"""
function UpdateTracers(Tracers, Grid, T, Phi, InterpolationMethod="Quadratic");

    dim = length(Grid)
    if isassigned(Tracers,1)        # only if the Tracers StructArray is non-empty
        # extract coordinates
        coord = Tracers.coord; coord = hcat(coord...)';       # extract array with coordinates of tracers

        x   = coord[:,1];
        z   = coord[:,end];
        if dim==2;
            Points_irregular = (x,z);
        else
            y   = coord[:,2];
            Points_irregular = (x,y,z);
        end

        # Correct coordinates (to stay withoin bounds of models)
        CorrectBounds!(Points_irregular, Grid);

        # Interpolate temperature from grid to tracers
        T_tracers           = tuple(zeros(size(x)));
        Interpolate!(T_tracers,         Grid, tuple(T), Points_irregular, InterpolationMethod);

        Phi_melt_tracers    = tuple(zeros(size(x)));
        Interpolate!(Phi_melt_tracers,  Grid, tuple(1.0 .- Phi), Points_irregular, InterpolationMethod);    # 1-Phi, as Phi=solid fraction

        # Update info on tracers
        for iT = 1:length(Tracers)
            LazyRow(Tracers, iT).T          = T_tracers[1][iT];             # Temperature
            LazyRow(Tracers, iT).Phi_melt   = Phi_melt_tracers[1][iT];      # Melt fraction
        end
    end

    return Tracers

end

"""
    This evenly populates the grid with tracers

        General form:

        Tracers     = InitializeTracers(Grid, NumTracersDir=3, RandomPertur=true);

        with:
            Grid:   2D or 3D arrays that describe the grid coordinates
                    2D - (X,Z)
                    3D - (X,Y,Z)

            NumTracersDir:    The number of tracers per direction

            RandomPertur:     Add slight random noise on tracer location?

        out:
            Tracers:    Tracers structure
"""
function InitializeTracers(Grid, NumTracersDir=3, RandomPertur=true);

    dim             =   length(Grid)
    R               =   CartesianIndices(Grid[1])
    Ifirst, Ilast   =   first(R), last(R)
    I1              =   oneunit(Ifirst)
    numTr           =   0;
    cen,d           =   zeros(dim), zeros(dim)
    coord_loc       =   zeros(NumTracersDir^dim,dim);
    Tracers         =   StructArray{Tracer}(undef, 1)                                    # Initialize Tracers structure

    # create a 'basic' Tracers struct
    t               =   Tracer(num=numTr, coord=zeros(dim), T=0.0, Phase=1);
    Tracers0        =   StructArray([t]);
    for i=1:NumTracersDir^dim-1
        append!(Tracers0,[t])
    end

    # The main assumption is that the grid coordinates specified in GRID are the corner points
    # of the cells. The control is given between the cells, as illustrated in the sketch below
    #
    #  X(ix,iz+1)
    #       | o  o  o  |
    #       x----------x
    #       |  o   o   |
    #       |   o    o |
    #  X(ix,iz)       X(ix+1,iz)
    #       x----------x
    #
    #
    # Note that the CartesianIndices in julia allow writing quite general code that works in any dimension

    # Determine center coordinate of current cell & width in every direction
    for idim=1:dim
        d[idim]    = (Grid[idim][Ifirst+I1] -   Grid[idim][Ifirst]);         # spacing of grid cells
    end
    d1      =   d/(NumTracersDir);
    xl      =   (-d[  1]/2. + d1[  1]/2.)
    zl      =   (-d[end]/2. + d1[end]/2.)
    xs      =   xl:d1[  1]: (xl+(NumTracersDir-1)*d1[  1]);
    zs      =   zl:d1[end]: (zl+(NumTracersDir-1)*d1[end]);
    if dim==3
        yl  =   (-d[2]/2. + d1[2]/2.)
        ys  =   yl:d1[2]: (yl+(NumTracersDir-1)*d1[2]);
    end

    # Generate local, regular, coordinate arrays for the new tracers
    if dim==2
        coord_loc0 = [ [x,z]     for x=xs for z=zs];             # Creates a tuple with coords
    elseif dim==3
        coord_loc0 = [ [x,y,z]   for x=xs for y=ys for z=zs];    # Creates a tuple with coords
    end
    coord_loc = coord_loc0;

    for I = Ifirst:Ilast-I1

        for idim=1:dim
            cen[idim]  = (Grid[idim][I+I1] +   Grid[idim][I])/2.0;     # center of control volume
            d[idim]    = (Grid[idim][I+I1] -   Grid[idim][I]);         # spacing of grid cells
        end

        ## THIS SECTION ALLOCATES A LOT EVEN WHEN WE HAVE IT FALSE  ---
        # add random perturbation if requested
        if RandomPertur

            if dim==2
                randm       = [(map(rand,(Float64,Float64)) .- 0.5).*2.0 .* (d1[1],d1[2]) for i=eachindex(coord_loc)];        #
                coord_loc   = [ coord_loc0[i] .+ randm[i] for i=eachindex(coord_loc0)];       # This can potentially be done in-place?
            elseif dim==3
                randm       = [(map(rand,(Float64,Float64, Float64)) .- 0.5).*2.0 .* (d1[1],d1[2],d1[3]) for i=eachindex(coord_loc)];        #
                coord_loc   = [ coord_loc0[i] .+ randm[i] for i=eachindex(coord_loc0)];       # This can potentially be done in-place?
            end
        end

        #if 1==0

        # Add new tracers with perturbed coords to struct

        ## THIS LINE ALLOCATES EVEN MORE, BECAUSE IT IS CALLED MANY TIMES
        ReplaceTracerFields!(Tracers0, coord_loc,cen,size(Tracers,1));     # Replace coord & num and add cen to coordinate

        ## THIS LINE IS SLOW BUT DOESN'T ALLOCATE ALL THAT MUCH:
        append!(Tracers, Tracers0);                                             # Extend the Tracers structure and add new fields to it
       # end
    end

    # delete first field, which was empty
    StructArrays.foreachfield(v -> deleteat!(v, 1), Tracers);

    return Tracers

end

function ReplaceTracerFields!(Tracers0, coord_loc::Array, cen::Array, num_start::Int)
    for i=1:length(coord_loc)
        Tracers0.coord[i]  = coord_loc[i] + cen;
        Tracers0.num[i]    = i+num_start-1;
    end
end

"""
    This computes the PhaseRatio on the grid specified by Grid

        General form:

        PhaseRatio,NumTracers = PhaseRatioFromTracers(Grid, Tracers, InterpolationMethod, BackgroundPhase=2, RequestNumTracers=false; BackgroundPhase=2);

        with:
            Grid:   2D or 3D arrays that describe the grid coordinates
                    2D - (X,Z)
                    3D - (X,Y,Z)

            Tracers:    Tracers structure

            InterpolationMethod:    Interpolation method used to go from Tracers ->  Grid
                    "Constant"          -   All particles within a distance [dx,dy,dz] around the grid point
                    "DistanceWeighted"  -   Particles closer to the grid point have a stronger weight.
                                            This follows what is described in:
                                                Duretz, T., May, D.A., Gerya, T.V., Tackley, P.J., 2011. Discretization errors and
                                                free surface stabilization in the finite difference and marker-in-cell method for applied geodynamics:
                                                A numerical study: Geochem. Geophys. Geosyst. 12, https://doi.org/10.1029/2011GC00356

            BackgroundPhase:        The background phase (used for places that don't have cells, nor surrounding cells )

            RequestNumTracers:      Return the number of tracers on every grid cell (default=false)

        optional parameters:

            BackgroundPhase:        if this is defined, points where no tracers are present, will get the number BackgroundPhase.
                                    You need to define this with a keywords as: (, BackgroundPhase=2)

        out:
            PhaseRatio:    Phase ratio on the gridpoints defined by Grid

            NumTracers:    The number of tracers per grid point

"""
function PhaseRatioFromTracers(FullGrid, Grid, Tracers, InterpolationMethod="Constant", RequestNumTracers=false; BackgroundPhase=nothing);

    numPhases       =   maximum(Tracers.Phase);
    if ~isnothing(BackgroundPhase)
        if numPhases<BackgroundPhase; numPhases=BackgroundPhase; end
    end
    dim             =   length(FullGrid)
    Nx              =   size(FullGrid[1],1);
    if dim==2
        Nz          =   size(FullGrid[1],2);
        siz         =   [Nx,Nz,numPhases];
    else
        Ny, Nz      =   size(FullGrid[1],2), size(FullGrid[1],3);
        siz         =   [Nx,Ny,Nz,numPhases];
    end
    PhaseRatio      =   zeros(Tuple(siz));                  # will hold phase ratio @ the end for every phase

    R               =   CartesianIndices(FullGrid[1])
    Ifirst, Ilast   =   first(R), last(R)
    I1              =   oneunit(Ifirst)

    # We assume that spacing is constant in all directions;
    #   If that is not the case the algorithm becomes a bit more complicated (not implemented)
    d               =   zeros(dim)
    for idim=1:dim
        d[idim]     =   (FullGrid[idim][Ifirst+I1] -   FullGrid[idim][Ifirst]);         # spacing of grid cells
    end
    coord           =   Tracers.coord; coord = hcat(coord...)';                         # extract array with coordinates of all tracers

    # Correct coordinates of tracers (to stay within bounds of grid), to not mess up the interpolation below
    CorrectBounds_Array!(coord, Grid);

    iPhase          =   Tracers.Phase;
    indX            =   CartesianIndices( FullGrid[  1]);

    IndPoints       =   zeros(Int,length(Tracers),1);
    if dim==2
        itp         =   interpolate(LinearIndices(FullGrid[1]), BSpline(Constant()));   # 2D interpolation that has indices
        sitp        =   scale(itp, Grid[1], Grid[2]);
        evaluate_interp_Int_2D(IndPoints, sitp,coord);                                  # this gives the full 2D/3D index
    else
        itp         =   interpolate(LinearIndices(FullGrid[1]), BSpline(Constant()));   # 2D interpolation that has indices
        sitp        =   scale(itp, Grid[1], Grid[2], Grid[3]);
        evaluate_interp_Int_3D(IndPoints, sitp,coord);                                  # this gives the full 2D/3D index
    end

    NumTracers  =   zeros(Int64,Tuple(siz[1:end-1]));   # Keep track of # of tracers around every point
    indNum      =   CartesianIndices(NumTracers);

    if InterpolationMethod=="DistanceWeighted"
        # In case we use a distance based weighting, compute the weight factor here
        X_g             =   FullGrid[1  ][IndPoints];  # Gridpoint to which the particle belongs
        Z_g             =   FullGrid[end][IndPoints];  # Gridpoint to which the particle belongs
        Weight      =   zeros(size(coord,1),1)
        if dim==2
            evaluate_weight_2D(Weight, coord, X_g, Z_g, d[1], d[end]);
        elseif dim==3
            Y_g         =   FullGrid[2  ][IndPoints];
            evaluate_weight_3D(Weight, coord, X_g, Y_g, Z_g, d[1], d[2], d[end]);
        end

    elseif InterpolationMethod=="Constant"
        Weight          =   ones(length(Tracers),1);
    else
        error("Unknown InterpolationMethod=$InterpolationMethod. Choose: [Constant] or [DistanceWeighted]. ")
    end

    indPhase             =   CartesianIndices(PhaseRatio);
    if dim==2
        Add_Phase_2D(PhaseRatio, NumTracers, length(Tracers), IndPoints, indPhase, indX, iPhase, indNum, Weight)
    elseif dim==3
        Add_Phase_3D(PhaseRatio, NumTracers, length(Tracers), IndPoints, indPhase, indX, iPhase, indNum, Weight)
    end

    if ~isnothing(BackgroundPhase)
        # All "empty" cells will be set to the background phase
        #
        # This is particularly useful in case we have models where Tracers do not 'fill'
        # the full model space

        ind_empty       = findall(x->x==0, NumTracers);
        for I in ind_empty
            # Use phase ratios from nearby point with most tracers
            if dim==2;  PhaseRatio[:,:  ,BackgroundPhase]    .=   1.0;
            else        PhaseRatio[:,:,:,BackgroundPhase]    .=   1.0;
            end

            NumTracers[I] = 1;
        end

    else
        # This is for cases where we do cover the full model domain
        #
        # Deal with points that have zero tracers @ this stage
        #  We do this by a query of the surrounding (9/27 in 2D/3D) points and use the properties of the closest point
        #
        # TODO: Deal with cases in which tracers are not defined in the full domain, but only used to track dikes
        #       and their properties, while otherwise having a fixed background
        ind_empty       = findall(x->x==0, NumTracers);
        R               = CartesianIndices(NumTracers)
        Ifirst, Ilast   = first(R), last(R)
        I1              = oneunit(Ifirst)
        for I in ind_empty
            @show I, NumTracers[I]

            # query surrounding points & store the index of the point with the most particles
            smax = 0;   Imax = I;
            for J in max(Ifirst, I-I1):min(Ilast, I+I1)
                if NumTracers[J]>smax
                    Imax = J;
                    smax = NumTracers[J];
                end
            end

            # Use phase ratios from nearby point with most tracers
            for iPhase=1:numPhases
                if dim==2;  PhaseRatio[I[1],I[2],     iPhase]   =   PhaseRatio[Imax[1],Imax[2],iPhase];
                else        PhaseRatio[I[1],I[2],I[3],iPhase]   =   PhaseRatio[Imax[1],Imax[2],Imax[3],iPhase];
                end
            end

            NumTracers[I] = 1;
        end

      # in case we still have problematic cells, at least spit out a warning
      if any(NumTracers .== 0 )
        warning("We still have cells that have no tracers and we did not manage to fix this by using nearby cells. ")
        end

    end
    sumPhaseRatio   =   sum(PhaseRatio,dims=dim+1);

    # Normalize phase ratio, such that sum=1
    for iPhase=1:numPhases
        if dim==2;  PhaseRatio[:,:  ,iPhase]    =   PhaseRatio[:,:  ,iPhase]./sumPhaseRatio;
        else        PhaseRatio[:,:,:,iPhase]    =   PhaseRatio[:,:,:,iPhase]./sumPhaseRatio;
        end
    end

    if RequestNumTracers
        return PhaseRatio, NumTracers
    else
        return PhaseRatio
    end
end

"""
    This averages a certain property from tracers -> grid. The data will ONLY be set on the grid if we have at least x tracers!
        if not, the origibal data will be kept

        General form:

            TracersToGrid!(Data, FullGrid, Grid, Tracers, Property="T", InterpolationMethod="Constant", RequestNumTracers=false);

        with:
            Data:       2D or 3D arrays that have the grid coordinates

            FullGrid:   2D or 3D arrays that describe the grid coordinates
                        2D - (X,Z)
                        3D - (X,Y,Z)

            Grid:       1D vectors describing the cartesian grid

            Tracers:    Tracers structure

            Property:   Property to be interpolated from Tracers2Grid
                        "T" - temperature

            InterpolationMethod:    Interpolation method used to go from Tracers ->  Grid
                    "Constant"          -   All particles within a distance [dx,dy,dz] around the grid point
                    "DistanceWeighted"  -   Particles closer to the grid point have a stronger weight.
                                            This follows what is described in:
                                                Duretz, T., May, D.A., Gerya, T.V., Tackley, P.J., 2011. Discretization errors and
                                                free surface stabilization in the finite difference and marker-in-cell method for applied geodynamics:
                                                A numerical study: Geochem. Geophys. Geosyst. 12, https://doi.org/10.1029/2011GC003567

            RequestNumTracers:    Return the number of tracers on every grid cell (default=false)

        out:
            PhaseRatio:    Phase ratio on the gridpoints defined by Grid

            NumTracers:    The number of tracers per grid point

"""
function TracersToGrid!(Data, FullGrid, Grid, Tracers, Property="T", InterpolationMethod="Constant", RequestNumTracers=false);

    numPhases       =   maximum(Tracers.Phase);
    dim             =   length(FullGrid)
    Nx              =   size(FullGrid[1],1);
    if dim==2
        Nz          =   size(FullGrid[1],2);
    else
        Ny, Nz      =   size(FullGrid[1],2), size(FullGrid[1],3);
    end
    Data_local      =   copy(Data) .* 0.0;                  # will hold the data @ the end

    R               =   CartesianIndices(FullGrid[1])
    Ifirst, Ilast   =   first(R), last(R)
    I1              =   oneunit(Ifirst)

    # We assume that spacing is constant in all directions;
    #   If that is not the case the algorithm becomes a bit more complicated (not implemented)
    d               =   zeros(dim)
    for idim=1:dim
        d[idim]     =   (FullGrid[idim][Ifirst+I1] -   FullGrid[idim][Ifirst]);         # spacing of grid cells
    end
    coord           =   Tracers.coord; coord = hcat(coord...)';                         # extract array with coordinates of all tracers

    # Correct coordinates of tracers (to stay within bounds of grid), to not mess up the interpolation below
    CorrectBounds_Array!(coord, Grid);

    indX            =   CartesianIndices( FullGrid[  1]);
    IndPoints       =   zeros(Int,length(Tracers),1);
    if dim==2
        itp         =   interpolate(LinearIndices(FullGrid[1]), BSpline(Constant()));   # 2D interpolation that has indices
        sitp        =   scale(itp, Grid[1], Grid[2]);
        evaluate_interp_Int_2D(IndPoints, sitp,coord);                                  # this gives the full 2D/3D index
    else
        itp         =   interpolate(LinearIndices(FullGrid[1]), BSpline(Constant()));   # 2D interpolation that has indices
        sitp        =   scale(itp, Grid[1], Grid[2], Grid[3]);
        evaluate_interp_Int_3D(IndPoints, sitp,coord);                                  # this gives the full 2D/3D index
    end

    NumTracers  =   zeros(Int64,  size(Data));   # Keep track of # of tracers around every point
    TotalWeight =   zeros(Float64,size(Data))
    indNum      =   CartesianIndices(NumTracers);

    if InterpolationMethod=="DistanceWeighted"
        # In case we use a distance based weighting, compute the weight factor here
        X_g             =   FullGrid[1  ][IndPoints];   # Gridpoint to which the particle belongs
        Z_g             =   FullGrid[end][IndPoints];   # Gridpoint to which the particle belongs
        Weight      =   zeros(size(coord,1),1)
        if dim==2
            evaluate_weight_2D(Weight, coord, X_g, Z_g, d[1], d[end]);
        elseif dim==3
            Y_g         =   FullGrid[2  ][IndPoints];
            evaluate_weight_3D(Weight, coord, X_g, Y_g, Z_g, d[1], d[2], d[end]);
        end

    elseif InterpolationMethod=="Constant"
        Weight          =   ones(length(Tracers),1);
    else
        error("Unknown InterpolationMethod=$InterpolationMethod. Choose: [Constant] or [DistanceWeighted]. ")
    end

    if  Property=="T"
        DataTracers          =   Tracers.T;     # temperature
    elseif Property=="Phi"
        DataTracers          =   Tracers.Phi;   # solid fraction
    else
        error("Property $(Property) not yet implemented")
    end

    SumData(Data_local, NumTracers, TotalWeight, length(Tracers), IndPoints, DataTracers, Weight);   # sum data from ever point

    # normalize based on weight & set non-empty data pounts
    Threads.@threads for i=eachindex(Data_local)
        if TotalWeight[i]>0.0
            Data_local[i] = Data_local[i]/TotalWeight[i];
        end
        if NumTracers[i]>0
            Data[i] = Data_local[i];
        end
    end

    if RequestNumTracers
        return NumTracers
    else
        return nothing
    end

end

# define functions to speed up key calculations above
function SumData(DataLocal, NumTracers, TotalWeight, nT, IndPoints, DataTracers, Weight)
    Threads.@threads    for iT=1:nT
                                ind = IndPoints[iT];
@inbounds                       DataLocal[ind]      +=  Weight[iT]*DataTracers[iT];
@inbounds                       TotalWeight[ind]    +=  Weight[iT];
@inbounds                       NumTracers[ind]     +=  1;
                        end
    end

function Add_Phase_2D(PhaseRatio, NumTracers, nT, IndPoints, indPhase, indX, iPhase, indNum, Weight)
    Threads.@threads    for iT=1:nT
                                iX = indX[IndPoints[iT]][1]
                                iZ = indX[IndPoints[iT]][2]
@inbounds                       PhaseRatio[   indPhase[iX,iZ,iPhase[iT]]]    +=  Weight[iT];
@inbounds                       NumTracers[indNum[iX,iZ             ]]       +=  1;
                        end
    end

function Add_Phase_3D(PhaseRatio, NumTracers, nT, IndPoints, indPhase, indX, iPhase, indNum, Weight)
        Threads.@threads    for iT=1:nT
                                    iX = indX[IndPoints[iT]][1]
                                    iY = indX[IndPoints[iT]][2]
                                    iZ = indX[IndPoints[iT]][3]
    @inbounds                       PhaseRatio[   indPhase[iX,iY,iZ,iPhase[iT]]]    +=  Weight[iT];
    @inbounds                       NumTracers[indNum[iX,iY,iZ             ]]       +=  1;
                            end
        end

function evaluate_weight_2D(w, coord, X_g, Z_g, dx, dz)
    Threads.@threads    for i=1:size(coord,1)
                           wx = 1.0 - abs.(coord[i,1  ] - X_g[i])/(dx/2.0);
                           wz = 1.0 - abs.(coord[i,2  ] - Z_g[i])/(dz/2.0);
        @inbounds          w[i]    = wx*wz;
                        end
    end

function evaluate_weight_3D(w, coord, X_g, Y_g, Z_g, dx, dy, dz)
        Threads.@threads    for i=1:size(coord,1)
                               wx = 1.0 - abs.(coord[i,1  ] - X_g[i])/(dx/2.0);
                               wy = 1.0 - abs.(coord[i,2  ] - Y_g[i])/(dy/2.0);
                               wz = 1.0 - abs.(coord[i,3  ] - Z_g[i])/(dz/2.0);
            @inbounds          w[i]    = wx*wy*wz;
                            end
        end
function evaluate_interp_Int_2D(s, itp, Points_irregular)
    Threads.@threads    for i=1:size(Points_irregular,1)
                            s[i]    = Int(itp(Points_irregular[i,1],Points_irregular[i,2]));
                        end
    end

function evaluate_interp_Int_3D(s, itp, Points_irregular)
    Threads.@threads    for i=1:size(Points_irregular,1)
                            s[i]    = Int(itp(Points_irregular[i,1],Points_irregular[i,2],Points_irregular[i,3]));
                        end
    end

"""
    CorrectBounds_Array!(Points, Grid);

Ensures that the coordinates of Points stay within the bounds
of the regular grid Grid. Points is an array of size [nPoints, dim]

 """
function CorrectBounds_Array!(Points, Grid);

    minC = [minimum(Grid[i]) for i=1:length(Grid)];
    maxC = [maximum(Grid[i]) for i=1:length(Grid)];

    Threads.@threads   for iT=1:size(Points,1)
        for i=1:length(Grid);
            if Points[iT,i] < minC[i]
                Points[iT, i ]  =   minC[i];
            end
            if Points[iT,i] > maxC[i]
                Points[iT, i ]  =   maxC[i];
            end
        end
    end
end

"""
        RockType = RockAssemblage(PhaseRatio)

        Computes the most abundant rock assemblage @ every point

"""
function RockAssemblage(PhaseRatio);
    # Most abundant rock-type @ every point
    dim         =   length(size(PhaseRatio))-1
    dummy       =   findmax(PhaseRatio,dims=dim+1);
    RockType    =   getindex.(dummy[2], dim+1)          # extract last of Cartesian index, which is the phase
    if dim==2
        RockType=RockType[:,:,1];
    end
    return RockType
end

"""
        CorrectTracersForTopography!(Tracers, Topo, PhaseAir=1, PhaseRock=2)

        Corrects tracers for topography, such that 'rock' tracers above the topography
            are set to "air" and vice-versa

        Input:
            Tracers     -       Tracers structure

            Topo        -
                                (Topo_x, Topo_y, Topo_z)  : topography in 3D
                                (Topo_x, Topo_z)  :         topography in 2D
            PhaseAir    -       Phase of air
            PhaseRock   -       Phase of rock

"""
function CorrectTracersForTopography!(Tracers, Topo, PhaseAir=1, PhaseRock=2)
    # Correct Tracers for topography:
    #       all Tracers with rock phase>1 are set to "air",
    #       all Tracesr with phase==1 is set to phase=2 ("rock")

    dim = length(Topo)
    if dim==2
        interp_linear = LinearInterpolation((Topo[1]), Topo[end],  extrapolation_bc = Line());
    elseif dim==3
        interp_linear = LinearInterpolation((Topo[1], Topo[2]), Topo[end],  extrapolation_bc = Line());
    end

    minT,maxT     = minimum(Topo[end]), maximum(Topo[end])

    for iT=1:length(Tracers)
        if      (Tracers[iT].coord[end]>maxT) & (Tracers.Phase[iT]>PhaseAir)
            Tracers.Phase[iT] = PhaseAir;               # air
        elseif (Tracers[iT].coord[end]<minT)  & (Tracers.Phase[iT]==PhaseAir)
            Tracers.Phase[iT] = PhaseRock;               # rock
        else
            # need to do interpolation
            pt      = Tracers[iT].coord;

            if dim==2
                z_topo  = interp_linear(pt[1]);
            elseif dim==3
                z_topo  = interp_linear(pt[1],pt[2]);
            end

            if      (pt[end]>z_topo) & (Tracers[iT].Phase>PhaseAir)
                Tracers.Phase[iT] = PhaseAir;          #
            elseif  (pt[end]<z_topo) & (Tracers[iT].Phase==PhaseAir)
                Tracers.Phase[iT] = PhaseRock;
            end

        end
    end

end

"""
        AdvectTracers!(Tracers, Grid, Velocity, dt, Method="RK2");

        Advects [Tracers] for one timestep (dt) using the [Velocity] defined on the points [Grid].

        Method: can be "Euler","RK2" or "RK4", for 1th, 2nd or 4th order explicit advection scheme, respectively.
"""
function AdvectTracers!(Tracers, Grid, Velocity, dt, Method="RK2");
    # Advect tracers forward in time & interpolate T on them

    dim             =   length(Grid);
    coord     =   hcat(Tracers.coord...)';    # extract array with coordinates of tracers

    x   = coord[:,1];
    z   = coord[:,end];
    if dim==2
        Points_irregular    =   (x,z);
    else
        y                   =   coord[:,2];
        Points_irregular    =   (x,y,z);
    end

    # Correct coordinates (to stay withoin bounds of models)
    CorrectBounds!(Points_irregular, Grid);

    # Advect
    Points_new              =   AdvectPoints(Points_irregular,  Grid,Velocity, dt,Method,  "Linear");     # Advect tracers

    # function to assign properties
    function testnoalloc_2D(sarr, val)
        for (Tracer,x,z) in zip(LazyRows(sarr), val[1], val[2])
            Tracer.coord = [x;z];
        end
    end

    function testnoalloc_3D(sarr, val)
        for (Tracer,x,y,z) in zip(LazyRows(sarr), val[1], val[2], val[3])
            Tracer.coord = [x; y; z];
        end
    end

    if dim==2
        testnoalloc_2D(Tracers, Points_new);
    else
        testnoalloc_3D(Tracers, Points_new);
    end

end


text:
#

using LinearAlgebra
using ProximalOperators
using ProximalAlgorithms.IterationTools
using Printf
using Base.Iterators
using Random
# using BregmanBC
using Flux

export solution

include("SGD_prox.jl")
include("SGD_prox_DNN.jl")

struct SGD{R<:Real}
    γ::Maybe{R}
    maxit::Int
    verbose::Bool
    freq::Int
    plus::Bool # true for diminishing stepsize
    DNN::Bool
    η0::Maybe{R}
    η_tilde::Maybe{R}
    function SGD{R}(;
        γ::Maybe{R} = nothing,
        maxit::Int = 10000,
        verbose::Bool = false,
        freq::Int = 1000,
        plus::Bool = false,
        DNN::Bool = false,
        η0::Maybe{R} = 0.1,
        η_tilde::Maybe{R} = 0.5,
    ) where {R}
        @assert γ === nothing || γ > 0
        @assert maxit > 0
        @assert freq > 0
        @assert η0 > 0
        @assert η_tilde > 0
        new(γ, maxit, verbose, freq, plus, DNN, η0, η_tilde)
    end
end

function (solver::SGD{R})(
    x0::AbstractArray{C};
    F = nothing,
    g = ProximalOperators.Zero(),
    L = nothing,
    μ = nothing,
    N = N,
) where {R,C<:RealOrComplex{R}}

    stop(state::SVRG_basic_state) = false
    disp(it, state) = @printf "%5d | %.3e  \n" it state.γ

    F === nothing && (F = fill(ProximalOperators.Zero(), (N,)))

    maxit = solver.maxit

    # dispatching the structure
    iter = SGD_prox_iterable(F, g, x0, N, L, μ, solver.γ)
    iter = take(halt(iter, stop), maxit)
    iter = enumerate(iter)
    num_iters, state_final = nothing, nothing
    for (it_, state_) in iter  # unrolling the iterator
        # see https://docs.julialang.org/en/v1/manual/interfaces/index.html
        if solver.verbose && mod(it_, solver.freq) == 0
            disp(it_, state_)
        end
        num_iters, state_final = it_, state_
    end
    if solver.verbose && mod(num_iters, solver.freq) !== 0
        disp(num_iters, state_final)
    end # for the final iteration
    return solution(state_final), num_iters
end

"""
    SGD([γ, maxit, verbose, freq])

Instantiate the SGD algorithm  for solving (strongly) convex optimization problems of the form

    minimize 1/N sum_{i=1}^N f_i(x) + g(x)

If `solver = SGD(args...)`, then the above problem is solved with

	solver(x0, [F, g, N, L, μ])

where F is an array containing f_i's, x0 is the initial point, and L, μ are arrays of
smoothness and strong convexity moduli of f_i's; they are optional when γ is provided.

Optional keyword arguments are:
* `γ`: stepsize
* `L`: an array of smoothness moduli of f_i's
* `μ`: (if strongly convex) an array of strong convexity moduli of f_i's
* `maxit::Integer` (default: `10000`), maximum number of iterations to perform.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `100`), frequency of verbosity.

"""

SGD(::Type{R}; kwargs...) where {R} = SGD{R}(; kwargs...)
SGD(; kwargs...) = SGD(Float64; kwargs...)

"""
If `solver = SVRG(args...)`, then

    itr = iterator(solver, x0, [F, g, N, L, μ])

is an iterable object. Note that [maxit, verbose, freq] fields of the solver are ignored here.

The solution at any given state can be obtained using solution(state), e.g.,
for state in Iterators.take(itr, maxit)
    # do something using solution(state)
end

See https://docs.julialang.org/en/v1/manual/interfaces/index.html
and https://docs.julialang.org/en/v1/base/iterators/ for a list of iteration utilities
"""

function iterator(
    solver::SGD{R},
    x0::Union{AbstractArray{C},Tp};
    F = nothing,
    g = ProximalOperators.Zero(),
    L = nothing,
    μ = nothing,
    N = N,
    data = nothing,
    DNN_config::Tdnn
) where {R,C<:RealOrComplex{R},Tp,Tdnn}
    F === nothing && (F = fill(ProximalOperators.Zero(), (N,)))
    # dispatching the iterator
    if solver.DNN
        iter = SGD_prox_DNN_iterable(F, g, x0, N, L, μ, solver.γ, solver.plus, solver.η0, solver.η_tilde, data, DNN_config)
    else
        iter = SGD_prox_iterable(F, g, x0, N, L, μ, solver.γ, solver.plus, solver.η0, solver.η_tilde)
    end
    return iter
end


text:
module IterationTools

export halt, tee, sample, stopwatch, loop

using Base.Iterators

# Halting

struct HaltingIterable{I, F}
    iter::I
    fun::F
end

Base.IteratorSize(::Type{HaltingIterable{I, F}}) where {I, F} = Base.IteratorSize(I)
Base.IteratorEltype(::Type{HaltingIterable{I, F}}) where {I, F} = Base.IteratorEltype(I)

Base.length(iter::HaltingIterable{I, F}) where {I, F} = length(iter.iter)
Base.eltype(iter::HaltingIterable{I, F}) where {I, F} = eltype(iter.iter)

function Base.iterate(iter::HaltingIterable)
    next = iterate(iter.iter)
    return dispatch(iter, next)
end

function Base.iterate(iter::HaltingIterable, (instruction, state))
    if instruction == :halt return nothing end
    next = iterate(iter.iter, state)
    return dispatch(iter, next)
end

function dispatch(iter::HaltingIterable, next)
    if next === nothing return nothing end
    return next[1], (iter.fun(next[1]) ? :halt : :continue, next[2])
end

halt(iter::I, fun::F) where {I, F} = HaltingIterable{I, F}(iter, fun)

# Side effects

struct TeeIterable{I, F}
    iter::I
    fun::F
end

Base.IteratorSize(::Type{TeeIterable{I, F}}) where {I, F} = Base.IteratorSize(I)
Base.IteratorEltype(::Type{TeeIterable{I, F}}) where {I, F} = Base.IteratorEltype(I)

Base.length(iter::TeeIterable{I, F}) where {I, F} = length(iter.iter)
Base.axes(iter::TeeIterable{I, F}) where {I, F} = axes(iter.iter)
Base.eltype(iter::TeeIterable{I, F}) where {I, F} = eltype(iter.iter)

function Base.iterate(iter::TeeIterable, args...)
    next = iterate(iter.iter, args...)
    if next !== nothing iter.fun(next[1]) end
    return next
end

tee(iter::I, fun::F) where {I, F} = TeeIterable{I, F}(iter, fun)

# Sampling

struct SamplingIterable{I}
    iter::I
    period::UInt
end

Base.IteratorSize(::Type{SamplingIterable{I}}) where I = Base.IteratorSize(I)
Base.IteratorEltype(::Type{SamplingIterable{I}}) where I = Base.IteratorEltype(I)

function Base.length(iter::SamplingIterable{I}) where I
    remainder = length(iter.iter) % iter.period
    quotient = length(iter.iter) ÷ iter.period
    return remainder == 0 ? quotient : quotient + 1
end

Base.size(iter::SamplingIterable{I}) where I = (length(iter), )
Base.eltype(iter::SamplingIterable{I}) where I = eltype(iter.iter)

function Base.iterate(iter::SamplingIterable, state=iter.iter)
    current = iterate(state)
    if current === nothing return nothing end
    for i = 1:iter.period-1
        next = iterate(state, current[2])
        if next === nothing return current[1], rest(state, current[2]) end
        current = next
    end
    return current[1], rest(state, current[2])
end

sample(iter::I, period) where I = SamplingIterable{I}(iter, period)

# Timing

struct StopwatchIterable{I}
    iter::I
end

Base.IteratorSize(::Type{StopwatchIterable{I}}) where I = Base.IteratorSize(I)
Base.IteratorEltype(::Type{StopwatchIterable{I}}) where I = Base.IteratorEltype(I)

Base.length(iter::StopwatchIterable{I}) where I = length(iter.iter)
Base.axes(iter::StopwatchIterable{I}) where I = axes(iter.iter)
Base.eltype(iter::StopwatchIterable{I}) where I = (UInt64, eltype(iter.iter))

function Base.iterate(iter::StopwatchIterable)
    t0 = time_ns()
    next = iterate(iter.iter)
    return dispatch(iter, t0, next)
end

function Base.iterate(iter::StopwatchIterable, (t0, state))
    next = iterate(iter.iter, state)
    return dispatch(iter, t0, next)
end

function dispatch(iter::StopwatchIterable, t0, next)
    if next === nothing return nothing end
    return (time_ns()-t0, next[1]), (t0, next[2])
end

stopwatch(iter::I) where I = StopwatchIterable{I}(iter)

# Looping

function loop(iter)
    x = nothing
    for y in iter x = y end
    return x
end

end # module


text:
import Winston
import Base.mimewritable

function save_winston(report, data)
  chunk = report.cur_chunk
  full_name, rel_name = get_figname(report, chunk)

  docformat = formats[report.formatdict[:doctype]]
  push!(report.figures, rel_name)
  report.fignum += 1
  vector_fmts = [".pdf"; ".svg"]
  #Don't use dpi for vector formats
  if chunk.options[:fig_ext] in vector_fmts
      Winston.savefig(data, full_name, width=chunk.options[:fig_width]*100,
          height=chunk.options[:fig_height]*100)
  else
      Winston.savefig(data, full_name,
          width=chunk.options[:fig_width]*chunk.options[:dpi],
          height=chunk.options[:fig_height]*chunk.options[:dpi])
  end
end

function Base.display(report::Report, m::MIME"image/svg+xml", data::Winston.FramedPlot)
  save_winston(report, data)
end

function Base.mimewritable(m::MIME"application/pdf", data::Winston.FramedPlot)
  return true
end

function Base.mimewritable(m::MIME"image/svg+xml", data::Winston.FramedPlot)
  return true
end

function Base.display(report::Report, m::MIME"application/pdf", data::Winston.FramedPlot)
  save_winston(report, data)
end

function Base.display(report::Report, m::MIME"image/png", data::Winston.FramedPlot)
  save_winston(report, data)
end

function Base.display(report::Report, m::MIME"image/jpg", data::Winston.FramedPlot)
  save_winston(report, data)
end


text:
function GaussLobattoQuad(n)
w=zeros(n+1,1);
x=zeros(n+1,1);
if n == 1
      w[1]=1;
      w[2]=1;
      x[1]=-1;
      x[2]=1;

elseif n == 2
      w[2]=4/3;
      w[1]=1/3;
      w[3]=1/3;

      x[2]=0;
      x[1]=-1;
      x[3]=1;

elseif n == 3
      w[2]=5/6;
      w[3]=5/6;
      w[1]=1/6;
      w[4]=1/6;

      x[2]=-1/sqrt(5);
      x[3]=1/sqrt(5);
      x[1]=-1;
      x[4]=1;

elseif n == 4
      w[3]=32/45;
      w[2]=49/90;
      w[4]=49/90;
      w[1]=1/10;
      w[5]=1/10;

      x[3]=0;
      x[2]=-sqrt(3/7);
      x[4]=sqrt(3/7);
      x[1]=-1;
      x[5]=1;

elseif n == 5
      w[3]=(14+sqrt(7))/30;
      w[4]=(14+sqrt(7))/30;
      w[2]=(14-sqrt(7))/30;
      w[5]=(14-sqrt(7))/30;
      w[1]=1/15;
      w[6]=1/15;

      x[3]=-sqrt(1/3-2*sqrt(7)/21);
      x[4]=sqrt(1/3-2*sqrt(7)/21);
      x[2]=-sqrt(1/3+2*sqrt(7)/21);
      x[5]=sqrt(1/3+2*sqrt(7)/21);
      x[1]=-1;
      x[6]=1;

elseif n == 6
      w[4]=256/525;
      w[3]=(124+7*sqrt(15))/350;
      w[5]=(124+7*sqrt(15))/350;
      w[2]=(124-7*sqrt(15))/350;
      w[6]=(124-7*sqrt(15))/350;
      w[1]=1/21;
      w[7]=1/21;

      x[4]=0;
      x[3]=-sqrt(5/11-2/11*sqrt(5/3));
      x[5]=sqrt(5/11-2/11*sqrt(5/3));
      x[2]=-sqrt(5/11+2/11*sqrt(5/3));
      x[6]=sqrt(5/11+2/11*sqrt(5/3));
      x[1]=-1;
      x[7]=1;

   else
      error("ord1+ord2 zu gro oder zu klein")
end
return (w,x)
end


text:
function ∇L_ρ_forward(f, gp::AbstractLatent, X::AbstractVector)
    θ, re = destructure(kernel(gp))
    return g = ForwardDiff.gradient(θ) do x
        k = re(x)
        Knn = kernelmatrix(k, X)
        f(Knn)
    end
end

function ∇L_ρ_forward(f, gp::SparseVarLatent, X::AbstractVector, ∇E_μ, ∇E_Σ, i, opt)
    θ, re = destructure(kernel(gp))
    return g = ForwardDiff.gradient(θ) do x
        k = re(x)
        Kmm = kernelmatrix(k, Zview(gp))
        Knm = kernelmatrix(k, X, Zview(gp))
        Knn = kernelmatrix_diag(k, X)
        f(Kmm, Knm, Knn, ∇E_μ, ∇E_Σ, i, opt)
    end
end

function ∇L_ρ_forward(f, gp::OnlineVarLatent, X::AbstractVector, ∇E_μ, ∇E_Σ, i, opt)
    θ, re = destructure(kernel(gp))
    return g = ForwardDiff.gradient(θ) do x
        k = re(x)
        Kmm = kernelmatrix(k, Zview(gp))
        Knm = kernelmatrix(k, X, Zview(gp))
        Knn = kernelmatrix_diag(k, X)
        Kaa = kernelmatrix(k, gp.Zₐ)
        Kab = kernelmatrix(k, gp.Zₐ, Zview(gp))
        f(Kmm, Knm, Knn, Kab, Kaa, ∇E_μ, ∇E_Σ, i, opt)
    end
end

## Return a function computing the gradient of the ELBO given the inducing point locations ##
function Z_gradient_forward(
    gp::SparseVarLatent{T},
    f_Z::Function,
    X,
    ∇E_μ::AbstractVector{T},
    ∇E_Σ::AbstractVector{T},
    i::AbstractInference,
    opt::InferenceOptimizer,
) where {T<:Real}
    gradient_inducing_points = similar(gp.Z.Z)
    #pre-allocation
    Jmm, Jnm = indpoint_derivative(kernel(gp), gp.Z),
    indpoint_derivative(kernel(gp), X, gp.Z)
    for j in 1:(gp.dim) #Iterate over the points
        for k in 1:size(gp.Z, 2) #iterate over the dimensions
            @views gradient_inducing_points[j, k] .= f_Z(
                Jmm[:, :, j, k], Jnm[:, :, j, k], ∇E_μ, ∇E_Σ, i, opt
            )
        end
    end
end

function Z_gradient_forward(
    gp::OnlineVarLatent{T},
    f_Z::Function,
    X,
    ∇E_μ::AbstractVector{T},
    ∇E_Σ::AbstractVector{T},
    i::AbstractInference,
    opt::InferenceOptimizer,
) where {T<:Real}
    Z_gradient = similar(gp.Z.Z)
    Jnm, Jab, Jmm = indpoint_derivative(kernel(gp), X, gp.Z),
    indpoint_derivative(kernel(gp), gp.Zₐ, gp.Z),
    indpoint_derivative(kernel(gp), gp.Z)
    for j in 1:(gp.dim) #Iterate over the points
        for k in 1:size(gp.Z, 2) #iterate over the dimensions
            @views Z_gradient[j, k] .= f_Z(
                Jmm[:, :, j, k], Jnm[:, :, j, k], Jab[:, :, j, k], ∇E_μ, ∇E_Σ, i, opt
            )
        end
    end
    return Z_gradient
end

function indpoint_derivative(kernel::Kernel, Z::AbstractVector)
    return reshape(
        ForwardDiff.jacobian(x -> kernelmatrix(kernel, x; obsdim=1), Z),
        size(Z, 1),
        size(Z, 1),
        size(Z, 1),
        size(Z, 2),
    )
end

function indpoint_derivative(kernel::Kernel, X, Z::AbstractVector)
    return reshape(
        ForwardDiff.jacobian(x -> kernelmatrix(kernel, X, x; obsdim=1), Z),
        size(X, 1),
        size(Z, 1),
        size(Z, 1),
        size(Z, 2),
    )
end


text:
using SparseArrays

function sub2ind(array_shape, rows, cols)
    out = (cols.-1).*array_shape[1] .+ rows
    return convert.(Int64, out)
end

function compute_warping_operator(flow)
    # ref: https://github.com/HendrikMuenster/flexBox/blob/79abc7285703911cca2653434cca5bcefc79c722/operators/generateMatrix/warpingOperator.m
    """
    u: gray image
    v: flow [H, W, 2]
    """
    H, W, dims = size(flow)
    shape = (H,W)
    XX, YY = repeat(collect(1:W)', H, 1), repeat(collect(1:H), 1, W)
    targetPoint2Mat = XX + flow[:,:,1]
    targetPoint1Mat = YY + flow[:,:,2]
    x1 = floor.(targetPoint1Mat).-1
    x2 = x1 .+ 1
    x3 = x1 .+ 2
    x4 = x1 .+ 3
    y1 = floor.(targetPoint2Mat).-1
    y2 = y1 .+ 1
    y3 = y1 .+ 2
    y4 = y1 .+ 3
    v2 = targetPoint1Mat .- x2
    v1 = targetPoint2Mat .- y2
    indicator = ones(H,W)
    indicator = (x1 .> 0).*indicator
    indicator = (x4 .<= W).*indicator
    indicator = (y1 .> 0).*indicator
    indicator = (y4 .<= H).*indicator
    indicator = indicator .> 0
    v1 = v1[indicator]
    v2 = v2[indicator]
    idxList = sub2ind(shape, YY, XX)
    idxList_ = collect(Iterators.flatten(idxList))
    #find the indices
    list_J  = vcat(
        idxList_[sub2ind(shape, x2[indicator], y2[indicator])],
        idxList_[sub2ind(shape, x3[indicator], y2[indicator])],
        idxList_[sub2ind(shape, x2[indicator], y3[indicator])],
        idxList_[sub2ind(shape, x3[indicator], y3[indicator])])
    # list_J  = vcat(
    #     idxList_[sub2ind(shape, x1[indicator], y1[indicator])],
    #     idxList_[sub2ind(shape, x2[indicator], y1[indicator])],
    #     idxList_[sub2ind(shape, x3[indicator], y1[indicator])],
    #     idxList_[sub2ind(shape, x4[indicator], y1[indicator])],
    #     idxList_[sub2ind(shape, x1[indicator], y2[indicator])],
    #     idxList_[sub2ind(shape, x2[indicator], y2[indicator])],
    #     idxList_[sub2ind(shape, x3[indicator], y2[indicator])],
    #     idxList_[sub2ind(shape, x4[indicator], y2[indicator])],
    #     idxList_[sub2ind(shape, x1[indicator], y3[indicator])],
    #     idxList_[sub2ind(shape, x2[indicator], y3[indicator])],
    #     idxList_[sub2ind(shape, x3[indicator], y3[indicator])],
    #     idxList_[sub2ind(shape, x4[indicator], y3[indicator])],
    #     idxList_[sub2ind(shape, x1[indicator], y4[indicator])],
    #     idxList_[sub2ind(shape, x2[indicator], y4[indicator])],
    #     idxList_[sub2ind(shape, x3[indicator], y4[indicator])],
    #     idxList_[sub2ind(shape, x4[indicator], y4[indicator])])
    list_J = collect(Iterators.flatten(list_J))
    list_val  = ( (1 .-v1) .* (1 .-v2))
    list_val = vcat(list_val,(v2 .* (1 .-v1)))
    list_val = vcat(list_val,(v1 .* (1 .-v2)))
    list_val = vcat(list_val,(v1 .* v2))
    # list_val =vcat(
    #     (v1.*v2.*(v1 .- 1).^2 .*(v2 .- 1).^2)/4,
    #     -1 .*(v1.*(v1 .- 1).^2 .*(3 .*v2.^3 .- 5 .*v2.^2 .+ 2))./4,
    #     -1 .*(v1.*v2.*(v1 .- 1).^2 .*(- 3 .*v2.^2 + 4 .*v2 .+ 1))./4,
    #     -1 .*(v1.*v2.^2 .*(v1 .- 1).^2 .*(v2 .- 1))./4,
    #     -1 .*(v2.*(v2 .- 1).^2 .*(3 .*v1.^3 .- 5 .*v1.^2 .+ 2))./4,
    #     ((3 .*v1.^3 .- 5 .*v1.^2 .+ 2).*(3 .*v2.^3 .- 5 .*v2.^2 .+ 2))./4,
    #     (v2.*(- 3 .*v2.^2 .+ 4 .*v2 .+ 1).*(3 .*v1.^3 .- 5 .*v1.^2 .+ 2))./4,
    #     (v2.^2 .*(v2 .- 1).*(3 .*v1.^3 .- 5 .*v1.^2 .+ 2))./4,
    #     -1 .*(v1.*v2.*(v2 .- 1).^2 .*(- 3 .*v1.^2 .+ 4 .*v1 .+ 1))./4,
    #     (v1.*(- 3 .*v1.^2 .+ 4 .*v1 .+ 1).*(3 .*v2.^3 .- 5 .*v2.^2 .+ 2))./4,
    #     (v1.*v2.*(- 3 .*v1.^2 .+ 4 .*v1 .+ 1).*(- 3 .*v2.^2 .+ 4 .*v2 .+ 1))./4,
    #     (v1.*v2.^2 .*(v2 .- 1).*(- 3 .*v1.^2 .+ 4 .*v1 .+ 1))./4,
    #     -1 .*(v1.^2 .*v2.*(v1 .- 1).*(v2 .- 1).^2)./4,
    #     (v1.^2 .*(v1 .- 1).*(3 .*v2.^3 .- 5 .*v2.^2 .+ 2))./4,
    #     (v1.^2 .*v2.*(v1 .- 1).*(- 3 .*v2.^2 .+ 4 .*v2 .+ 1))./4,
    #     (v1.^2 .*v2.^2 .*(v1 .- 1).*(v2 .- 1))./4
    # )
    list_val = collect(Iterators.flatten(list_val))
    list_I = idxList[indicator]
    #list_I =collect(Iterators.flatten(repeat(list_I,1,16)))
    list_I =collect(Iterators.flatten(repeat(list_I,1,4)))
    # list_I_py, list_J_py, list_val_py = py"py_warping_operator"(flow)
    # list_I_py = list_I_py.+1
    # list_J_py = list_J_py.+1
    # W_mat_py = sparse(list_I_py, list_J_py, list_val_py, H*W, H*W)
    W_mat = dropzeros!(sparse(list_I, list_J, list_val, H*W, H*W))
    return W_mat
end

img2 = zeros(128,128)
H, W = size(img2)
img2[30:60,40:50] .= 255.0

flow = zeros(H, W, 2)
flow[30:60,40:50, 1] .= -1.0

Wop = compute_warping_operator(flow)
    img1_ = Wop * vec(img2)
img1 = reshape(img1_, H, W)

max_stars_repo_name: QPImaging/xfromProjections

max_issues_repo_name: QPImaging/xfromProjections

max_forks_repo_name: QPImaging/xfromProjections


text:
using Documenter, ExcelReaders

makedocs(
	modules = [ExcelReaders],
	sitename = "ExcelReaders.jl",
	analytics="UA-132838790-1",
	pages = [
        "Introduction" => "index.md"
    ]
)

deploydocs(
    repo = "github.com/queryverse/ExcelReaders.jl.git"
)


text:
# ---
# title: 21. Merge Two Sorted Lists
# id: problem21
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/merge-two-sorted-lists/description/>
# hidden: true
# ---
#
# Merge two sorted linked lists and return it as a new **sorted** list. The new
# list should be made by splicing together the nodes of the first two lists.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/03/merge_ex1.jpg)
#
#
#
#     Input: l1 = [1,2,4], l2 = [1,3,4]
#     Output: [1,1,2,3,4,4]
#
#
# **Example 2:**
#
#
#
#     Input: l1 = [], l2 = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: l1 = [], l2 = [0]
#     Output: [0]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in both lists is in the range `[0, 50]`.
#   * `-100 <= Node.val <= 100`
#   * Both `l1` and `l2` are sorted in **non-decreasing** order.
#
#
## @lc code=start
using LeetCode

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


text:
function _laplace_diag_kernel!(
    dst     ::CuDeviceVector{T},
    Ξ       ::CuPositionVector{T},
    elements::CuTriangleVector{T},
    pot     ::F
) where {T, F <: Function}
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    i > length(Ξ) && return

    dst[i] = pot(Ξ[i], elements[i])
    nothing
end

for L in [:single, :double]
    @eval begin
        @inline $(Symbol("_laplace_", L, "_diag_kernel!"))(
            dst     ::CuDeviceVector{T},
            Ξ       ::CuPositionVector{T},
            elements::CuTriangleVector{T}
        ) where T = _laplace_diag_kernel!(dst, Ξ, elements, $(Symbol("laplacepot_", L)))
    end
end

function _laplace_mul_kernel!(
    dst     ::CuDeviceVector{T},
    Ξ       ::CuPositionVector{T},
    elements::CuTriangleVector{T},
    x       ::CuDeviceVector{T},
    pot     ::F
) where {T, F <: Function}
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    i > length(Ξ) && return

    ξ = Ξ[i]
    val = T(0)
    for j in 1:length(elements)
        val = CUDA.fma(pot(ξ, elements[j]), x[j], val)
    end
    dst[i] = val
    nothing
end

for L in [:single, :double]
    @eval begin
        @inline $(Symbol("_laplace_", L, "_mul_kernel!"))(
            dst     ::CuDeviceVector{T},
            Ξ       ::CuPositionVector{T},
            elements::CuTriangleVector{T},
            x       ::CuDeviceVector{T}
        ) where T = _laplace_mul_kernel!(dst, Ξ, elements, x, $(Symbol("laplacepot_", L)))
    end
end

function _reyukawa_diag_kernel!(
    dst     ::CuDeviceVector{T},
    Ξ       ::CuPositionVector{T},
    elements::CuTriangleVector{T},
    yuk     ::T,
    pot     ::F
) where {T, F <: Function}
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    i > length(Ξ) && return

    dst[i] = pot(Ξ[i], elements[i], yuk)
    nothing
end

for L in [:single, :double]
    @eval begin
        @inline $(Symbol("_reyukawa_", L, "_diag_kernel!"))(
            dst     ::CuDeviceVector{T},
            Ξ       ::CuPositionVector{T},
            elements::CuTriangleVector{T},
            yuk     ::T
        ) where T = _reyukawa_diag_kernel!(dst, Ξ, elements, yuk, $(Symbol("regularyukawapot_", L)))
    end
end

function _reyukawa_mul_kernel!(
    dst     ::CuDeviceVector{T},
    Ξ       ::CuPositionVector{T},
    elements::CuTriangleVector{T},
    x       ::CuDeviceVector{T},
    yuk     ::T,
    pot     ::F
) where {T, F <: Function}
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    i > length(Ξ) && return

    ξ = Ξ[i]
    val = T(0)
    for j in 1:length(elements)
        val = CUDA.fma(pot(ξ, elements[j], yuk), x[j], val)
    end
    dst[i] = val
    nothing
end

for L in [:single, :double]
    @eval begin
        @inline $(Symbol("_reyukawa_", L, "_mul_kernel!"))(
            dst     ::CuDeviceVector{T},
            Ξ       ::CuPositionVector{T},
            elements::CuTriangleVector{T},
            x       ::CuDeviceVector{T},
            yuk     ::T
        ) where T = _reyukawa_mul_kernel!(
            dst, Ξ, elements, x, yuk, $(Symbol("regularyukawapot_", L))
        )
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

import Base: copy, adjoint, getindex, show, transpose, one, zero, inv,
             hcat, vcat, hvcat

"""
    UniformScaling{T<:Number}

Generically sized uniform scaling operator defined as a scalar times the
identity operator, `λ*I`. See also [`I`](@ref).

# Examples
```jldoctest
julia> J = UniformScaling(2.)
UniformScaling{Float64}
2.0*I

julia> A = [1. 2.; 3. 4.]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> J*A
2×2 Array{Float64,2}:
 2.0  4.0
 6.0  8.0
```
"""
struct UniformScaling{T<:Number}
    λ::T
end

"""
    I

An object of type [`UniformScaling`](@ref), representing an identity matrix of any size.

# Examples
```jldoctest
julia> fill(1, (5,6)) * I == fill(1, (5,6))
true

julia> [1 2im 3; 1im 2 3] * I
2×3 Array{Complex{Int64},2}:
 1+0im  0+2im  3+0im
 0+1im  2+0im  3+0im
```
"""
const I = UniformScaling(true)

eltype(::Type{UniformScaling{T}}) where {T} = T
ndims(J::UniformScaling) = 2
Base.has_offset_axes(::UniformScaling) = false
getindex(J::UniformScaling, i::Integer,j::Integer) = ifelse(i==j,J.λ,zero(J.λ))

function show(io::IO, J::UniformScaling)
    s = "$(J.λ)"
    if occursin(r"\w+\s*[\+\-]\s*\w+", s)
        s = "($s)"
    end
    print(io, "$(typeof(J))\n$s*I")
end
copy(J::UniformScaling) = UniformScaling(J.λ)

conj(J::UniformScaling) = UniformScaling(conj(J.λ))

transpose(J::UniformScaling) = J
adjoint(J::UniformScaling) = UniformScaling(conj(J.λ))

one(::Type{UniformScaling{T}}) where {T} = UniformScaling(one(T))
one(J::UniformScaling{T}) where {T} = one(UniformScaling{T})
oneunit(::Type{UniformScaling{T}}) where {T} = UniformScaling(oneunit(T))
oneunit(J::UniformScaling{T}) where {T} = oneunit(UniformScaling{T})
zero(::Type{UniformScaling{T}}) where {T} = UniformScaling(zero(T))
zero(J::UniformScaling{T}) where {T} = zero(UniformScaling{T})

istriu(::UniformScaling) = true
istril(::UniformScaling) = true
issymmetric(::UniformScaling) = true
ishermitian(J::UniformScaling) = isreal(J.λ)

(+)(J::UniformScaling, x::Number) = J.λ + x
(+)(x::Number, J::UniformScaling) = x + J.λ
(-)(J::UniformScaling, x::Number) = J.λ - x
(-)(x::Number, J::UniformScaling) = x - J.λ

(+)(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ+J2.λ)
(+)(B::BitArray{2}, J::UniformScaling)      = Array(B) + J
(+)(J::UniformScaling, B::BitArray{2})      = J + Array(B)
(+)(J::UniformScaling, A::AbstractMatrix)   = A + J

(-)(J::UniformScaling)                      = UniformScaling(-J.λ)
(-)(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ-J2.λ)
(-)(B::BitArray{2}, J::UniformScaling)      = Array(B) - J
(-)(J::UniformScaling, B::BitArray{2})      = J - Array(B)
(-)(A::AbstractMatrix, J::UniformScaling)   = A + (-J)

# Unit{Lower/Upper}Triangular matrices become {Lower/Upper}Triangular under
# addition with a UniformScaling
for (t1, t2) in ((:UnitUpperTriangular, :UpperTriangular),
                 (:UnitLowerTriangular, :LowerTriangular))
    @eval begin
        function (+)(UL::$t1, J::UniformScaling)
            ULnew = copy_oftype(UL.data, Base._return_type(+, Tuple{eltype(UL), typeof(J)}))
            for i in axes(ULnew, 1)
                ULnew[i,i] = one(ULnew[i,i]) + J
            end
            return ($t2)(ULnew)
        end
    end
end

function (+)(A::AbstractMatrix, J::UniformScaling)
    checksquare(A)
    B = copy_oftype(A, Base._return_type(+, Tuple{eltype(A), typeof(J)}))
    @inbounds for i in axes(A, 1)
        B[i,i] += J
    end
    return B
end

function (-)(J::UniformScaling, A::AbstractMatrix)
    checksquare(A)
    B = convert(AbstractMatrix{Base._return_type(+, Tuple{eltype(A), typeof(J)})}, -A)
    @inbounds for i in axes(A, 1)
        B[i,i] += J
    end
    return B
end

inv(J::UniformScaling) = UniformScaling(inv(J.λ))
opnorm(J::UniformScaling, p::Real=2) = opnorm(J.λ, p)

function det(J::UniformScaling{T}) where T
    if isone(J.λ)
        one(T)
    elseif iszero(J.λ)
        zero(T)
    else
        throw(ArgumentError("Determinant of UniformScaling is only well-defined when λ = 0 or 1."))
    end
end

*(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ*J2.λ)
*(B::BitArray{2}, J::UniformScaling) = *(Array(B), J::UniformScaling)
*(J::UniformScaling, B::BitArray{2}) = *(J::UniformScaling, Array(B))
*(A::AbstractMatrix, J::UniformScaling) = A*J.λ
*(J::UniformScaling, A::AbstractVecOrMat) = J.λ*A
*(x::Number, J::UniformScaling) = UniformScaling(x*J.λ)
*(J::UniformScaling, x::Number) = UniformScaling(J.λ*x)

/(J1::UniformScaling, J2::UniformScaling) = J2.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ/J2.λ)
/(J::UniformScaling, A::AbstractMatrix) = lmul!(J.λ, inv(A))
/(A::AbstractMatrix, J::UniformScaling) = J.λ == 0 ? throw(SingularException(1)) : A/J.λ

/(J::UniformScaling, x::Number) = UniformScaling(J.λ/x)

\(J1::UniformScaling, J2::UniformScaling) = J1.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ\J2.λ)
\(A::Union{Bidiagonal{T},AbstractTriangular{T}}, J::UniformScaling) where {T<:Number} =
    rmul!(inv(A), J.λ)
\(J::UniformScaling, A::AbstractVecOrMat) = J.λ == 0 ? throw(SingularException(1)) : J.λ\A
\(A::AbstractMatrix, J::UniformScaling) = rmul!(inv(A), J.λ)

\(x::Number, J::UniformScaling) = UniformScaling(x\J.λ)

Broadcast.broadcasted(::typeof(*), x::Number,J::UniformScaling) = UniformScaling(x*J.λ)
Broadcast.broadcasted(::typeof(*), J::UniformScaling,x::Number) = UniformScaling(J.λ*x)

Broadcast.broadcasted(::typeof(/), J::UniformScaling,x::Number) = UniformScaling(J.λ/x)

==(J1::UniformScaling,J2::UniformScaling) = (J1.λ == J2.λ)

## equality comparison with UniformScaling
==(J::UniformScaling, A::AbstractMatrix) = A == J
function ==(A::AbstractMatrix, J::UniformScaling)
    @assert !has_offset_axes(A)
    size(A, 1) == size(A, 2) || return false
    iszero(J.λ) && return iszero(A)
    isone(J.λ) && return isone(A)
    return A == J.λ*one(A)
end
function ==(A::StridedMatrix, J::UniformScaling)
    size(A, 1) == size(A, 2) || return false
    iszero(J.λ) && return iszero(A)
    isone(J.λ) && return isone(A)
    for j in axes(A, 2), i in axes(A, 1)
        ifelse(i == j, A[i, j] == J.λ, iszero(A[i, j])) || return false
    end
    return true
end

function isapprox(J1::UniformScaling{T}, J2::UniformScaling{S};
            atol::Real=0, rtol::Real=Base.rtoldefault(T,S,atol), nans::Bool=false) where {T<:Number,S<:Number}
    isapprox(J1.λ, J2.λ, rtol=rtol, atol=atol, nans=nans)
end
function isapprox(J::UniformScaling, A::AbstractMatrix;
                  atol::Real = 0,
                  rtol::Real = Base.rtoldefault(promote_leaf_eltypes(A), eltype(J), atol),
                  nans::Bool = false, norm::Function = norm)
    n = checksquare(A)
    normJ = norm === opnorm             ? abs(J.λ) :
            norm === LinearAlgebra.norm ? abs(J.λ) * sqrt(n) :
                                          norm(Diagonal(fill(J.λ, n)))
    return norm(A - J) <= max(atol, rtol * max(norm(A), normJ))
end
isapprox(A::AbstractMatrix, J::UniformScaling; kwargs...) = isapprox(J, A; kwargs...)

function copyto!(A::AbstractMatrix, J::UniformScaling)
    @assert !has_offset_axes(A)
    size(A,1)==size(A,2) || throw(DimensionMismatch("a UniformScaling can only be copied to a square matrix"))
    fill!(A, 0)
    λ = J.λ
    for i = 1:size(A,1)
        @inbounds A[i,i] = λ
    end
    return A
end

function cond(J::UniformScaling{T}) where T
    onereal = inv(one(real(J.λ)))
    return J.λ ≠ zero(T) ? onereal : oftype(onereal, Inf)
end

# promote_to_arrays(n,k, T, A...) promotes any UniformScaling matrices
# in A to matrices of type T and sizes given by n[k:end].  n is an array
# so that the same promotion code can be used for hvcat.  We pass the type T
# so that we can re-use this code for sparse-matrix hcat etcetera.
promote_to_arrays_(n::Int, ::Type{Matrix}, J::UniformScaling{T}) where {T} = copyto!(Matrix{T}(undef, n,n), J)
promote_to_arrays_(n::Int, ::Type, A::AbstractVecOrMat) = A
promote_to_arrays(n,k, ::Type) = ()
promote_to_arrays(n,k, ::Type{T}, A) where {T} = (promote_to_arrays_(n[k], T, A),)
promote_to_arrays(n,k, ::Type{T}, A, B) where {T} =
    (promote_to_arrays_(n[k], T, A), promote_to_arrays_(n[k+1], T, B))
promote_to_arrays(n,k, ::Type{T}, A, B, C) where {T} =
    (promote_to_arrays_(n[k], T, A), promote_to_arrays_(n[k+1], T, B), promote_to_arrays_(n[k+2], T, C))
promote_to_arrays(n,k, ::Type{T}, A, B, Cs...) where {T} =
    (promote_to_arrays_(n[k], T, A), promote_to_arrays_(n[k+1], T, B), promote_to_arrays(n,k+2, T, Cs...)...)
promote_to_array_type(A::Tuple{Vararg{Union{AbstractVecOrMat,UniformScaling}}}) = Matrix

for (f,dim,name) in ((:hcat,1,"rows"), (:vcat,2,"cols"))
    @eval begin
        function $f(A::Union{AbstractVecOrMat,UniformScaling}...)
            n = 0
            for a in A
                if !isa(a, UniformScaling)
                    @assert !has_offset_axes(a)
                    na = size(a,$dim)
                    n > 0 && n != na &&
                        throw(DimensionMismatch(string("number of ", $name,
                            " of each array must match (got ", n, " and ", na, ")")))
                    n = na
                end
            end
            n == 0 && throw(ArgumentError($("$f of only UniformScaling objects cannot determine the matrix size")))
            return $f(promote_to_arrays(fill(n,length(A)),1, promote_to_array_type(A), A...)...)
        end
    end
end

function hvcat(rows::Tuple{Vararg{Int}}, A::Union{AbstractVecOrMat,UniformScaling}...)
    @assert !has_offset_axes(A...)
    nr = length(rows)
    sum(rows) == length(A) || throw(ArgumentError("mismatch between row sizes and number of arguments"))
    n = zeros(Int, length(A))
    needcols = false # whether we also need to infer some sizes from the column count
    j = 0
    for i = 1:nr # infer UniformScaling sizes from row counts, if possible:
        ni = 0 # number of rows in this block-row
        for k = 1:rows[i]
            if !isa(A[j+k], UniformScaling)
                na = size(A[j+k], 1)
                ni > 0 && ni != na &&
                    throw(DimensionMismatch("mismatch in number of rows"))
                ni = na
            end
        end
        if ni > 0
            for k = 1:rows[i]
                n[j+k] = ni
            end
        else # row consisted only of UniformScaling objects
            needcols = true
        end
        j += rows[i]
    end
    if needcols # some sizes still unknown, try to infer from column count
        nc = j = 0
        for i = 1:nr
            nci = 0
            rows[i] > 0 && n[j+1] == 0 && continue # column count unknown in this row
            for k = 1:rows[i]
                nci += isa(A[j+k], UniformScaling) ? n[j+k] : size(A[j+k], 2)
            end
            nc > 0 && nc != nci && throw(DimensionMismatch("mismatch in number of columns"))
            nc = nci
            j += rows[i]
        end
        nc == 0 && throw(ArgumentError("sizes of UniformScalings could not be inferred"))
        j = 0
        for i = 1:nr
            if rows[i] > 0 && n[j+1] == 0 # this row consists entirely of UniformScalings
                nci = nc ÷ rows[i]
                nci * rows[i] != nc && throw(DimensionMismatch("indivisible UniformScaling sizes"))
                for k = 1:rows[i]
                    n[j+k] = nci
                end
            end
            j += rows[i]
        end
    end
    return hvcat(rows, promote_to_arrays(n,1, promote_to_array_type(A), A...)...)
end

## Matrix construction from UniformScaling
function Matrix{T}(s::UniformScaling, dims::Dims{2}) where {T}
    A = zeros(T, dims)
    v = T(s.λ)
    for i in diagind(dims...)
        @inbounds A[i] = v
    end
    return A
end
Matrix{T}(s::UniformScaling, m::Integer, n::Integer) where {T} = Matrix{T}(s, Dims((m, n)))
Matrix(s::UniformScaling, m::Integer, n::Integer) = Matrix(s, Dims((m, n)))
Matrix(s::UniformScaling, dims::Dims{2}) = Matrix{eltype(s)}(s, dims)
Array{T}(s::UniformScaling, dims::Dims{2}) where {T} = Matrix{T}(s, dims)
Array{T}(s::UniformScaling, m::Integer, n::Integer) where {T} = Matrix{T}(s, m, n)
Array(s::UniformScaling, m::Integer, n::Integer) = Matrix(s, m, n)
Array(s::UniformScaling, dims::Dims{2}) = Matrix(s, dims)

## Diagonal construction from UniformScaling
Diagonal{T}(s::UniformScaling, m::Integer) where {T} = Diagonal{T}(fill(T(s.λ), m))
Diagonal(s::UniformScaling, m::Integer) = Diagonal{eltype(s)}(s, m)

max_stars_repo_name: pochoi/julia

max_issues_repo_name: amilsted/julia

max_forks_repo_name: amilsted/julia


text:
using SlepcWrap
using Test

# from :
# https://discourse.julialang.org/t/what-general-purpose-commands-do-you-usually-end-up-adding-to-your-projects/4889
@generated function compare_struct(x, y)
    if !isempty(fieldnames(x)) && x == y
        mapreduce(n -> :(x.$n == y.$n), (a,b)->:($a && $b), fieldnames(x))
    else
        :(x == y)
    end
end

@testset "SlepcWrap.jl" begin
    include("./helmholtz.jl")
    include("./helmholtz_fancy.jl")
end