text:
"""
`module JAC.RadialIntegrals`
    ... a submodel of JAC that contains methods for calculating radial one- and two-particle matrix elements. These integrals occur
        frequently in atomic structure and collision theory, and their fast computations often appears essential.
"""
module  RadialIntegrals

    using  GSL, QuadGK
    using  ..AngularMomentum, ..Basics, ..Bsplines, ..Defaults,  ..Radial, ..Math, ..Nuclear
    ##x global JAC_counter = 0

    """
    `RadialIntegrals.GrantIab(a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid, potential::Radial.Potential)`
        ... computes the (radial) single-electron energy integral:

            `I(ab) = <a | h_D | b> = delta_{kappa_a, kappa_b} int_0^infty dr  [ c Q_a ( d/dr + kappa_a/r ) P_b +  c P_a (-d/dr + kappa_a/r ) Q_b
                                                                                - 2c^2 Q_a Q_b + V_nuc (r) (P_a P_b + Q_a Q_b) ]`

            for the orbitals a and b on the grid. potential.Zr must provide the effective nuclear charge Z(r) on this grid.
    """
    function GrantIab(a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid, potential::Radial.Potential)
        if  a.subshell.kappa != b.subshell.kappa    return( 0 )    end
        kappa = a.subshell.kappa;                   Zr = potential.Zr
        mtp   = min(size(a.P, 1), size(b.P, 1));    wc = Defaults.getDefaults("speed of light: c")

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp

            function f1(i::Int64)   dPb(j) = Math.derivative(b.P, j)
                                    dQb(j) = Math.derivative(b.Q, j)
                                    return( a.Q[i] * dPb(i) - a.P[i] * dQb(i) )                          end
            function f2(i::Int64)   return( (a.Q[i] * b.P[i] + a.P[i] * b.Q[i]) / grid.r[i] )            end
            function f3(i::Int64)   return( a.Q[i] * b.Q[i] )                                            end
            function f4(i::Int64)   return( -Zr[i] * (a.P[i] * b.P[i] + a.Q[i] * b.Q[i]) / grid.r[i] )   end

            I1 = Math.integrateFit(f1, min(size(a.P, 1), size(b.P, 1)), grid) / grid.h
            I2 = Math.integrateFitTransform(f2, min(size(a.P, 1), size(b.P, 1)), grid)
            I3 = Math.integrateFitTransform(f3, min(size(a.P, 1), size(b.P, 1)), grid)
            I4 = Math.integrateFitTransform(f4, min(size(a.P, 1), size(b.P, 1)), grid)

            return( Defaults.INVERSE_FINE_STRUCTURE_CONSTANT * I1 + Defaults.INVERSE_FINE_STRUCTURE_CONSTANT * kappa * I2 -
                    2 * Defaults.INVERSE_FINE_STRUCTURE_CONSTANT^2 * I3 + I4 )

        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp
               wa = wa + grid.wr[i] * (  wc * a.Q[i] * (b.Pprime[i] + kappa/grid.r[i] * b.P[i])
                                       - wc * a.P[i] * (b.Qprime[i] - kappa/grid.r[i] * b.Q[i])
                                       - 2wc^2 * a.Q[i] * b.Q[i]
                                       - Zr[i] * (a.P[i] * b.P[i] + a.Q[i] * b.Q[i]) / grid.r[i]  )
            end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.GrantILminus(L::Int64, q::Float64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)`
        ... computes Grant's (radial) integral for two relativistic orbitals:
            I_L^- (q; a,b) = int_0^\\infty dr j_L (qr) [ P_a Q_b - Q_a P_b ] .
    """
    function GrantILminus(L::Int64, q::Float64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)
        mtp = min(size(a.P, 1), size(b.P, 1))

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            function f(i :: Int64)    return( (a.P[i] * b.Q[i] - a.Q[i] * b.P[i]) * GSL.sf_bessel_jl(L, q * grid.r[i]) )       end
            return( Math.integrateFitTransform(f, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp   wa = wa + (a.P[i] * b.Q[i] - a.Q[i] * b.P[i]) * GSL.sf_bessel_jl(L, q * grid.r[i]) * grid.wr[i]   end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.GrantILplus(L::Int64, q::Float64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)`
        ... computes Grant's (radial) integral for two relativistic orbitals:
            I_L^+ (q; a,b) = int_0^\\infty dr j_L (qr) [ P_a Q_b + Q_a P_b ] .
    """
    function GrantILplus(L::Int64, q::Float64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)
        mtp = min(size(a.P, 1), size(b.P, 1))

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            function f(i :: Int64)    return( (a.P[i] * b.Q[i] + a.Q[i] * b.P[i]) * GSL.sf_bessel_jl(L, q * grid.r[i]) )       end
            return( Math.integrateFitTransform(f, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp   wa = wa + (a.P[i] * b.Q[i] + a.Q[i] * b.P[i]) * GSL.sf_bessel_jl(L, q * grid.r[i]) * grid.wr[i]   end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.GrantIL0(L::Int64, q::Float64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)`
        ... computes Grant's (radial) integral for two relativistic orbitals:
            I_L^0 (q; a,b) = int_0^\\infty dr j_L (qr) [ P_a Q_b ] .
    """
    function GrantIL0(L::Int64, q::Float64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)
        mtp = min(size(a.P, 1), size(b.P, 1))

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            function f(i :: Int64)    return( (a.P[i] * b.Q[i]) * GSL.sf_bessel_jl(L, q * grid.r[i]) )       end
            return( Math.integrateFitTransform(f, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp   wa = wa + (a.P[i] * b.Q[i]) * GSL.sf_bessel_jl(L, q * grid.r[i]) * grid.wr[i]   end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.GrantJL(L::Int64, q::Float64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)`
        ... computes Grant's (radial) integral for two relativistic orbitals:
            J_L (q; a,b) = int_0^\\infty dr j_L (qr) [ P_a P_b + Q_a Q_b ] .
    """
    function GrantJL(L::Int64, q::Float64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)
        mtp = min(size(a.P, 1), size(b.P, 1))

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            function f(i :: Int64)    return( (a.P[i] * b.P[i] + a.Q[i] * b.Q[i]) * GSL.sf_bessel_jl(L, q * grid.r[i]) )       end
            return( Math.integrateFitTransform(f, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp   wa = wa + (a.P[i] * b.P[i] + a.Q[i] * b.Q[i]) * GSL.sf_bessel_jl(L, q * grid.r[i]) * grid.wr[i]   end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.isotope_nms(a::Orbital, b::Orbital, Z::Float64, grid::Radial.Grid)`
        ... computes the normal mass shift radial integral int_o^infty ... A value::Float64 is returned.
    """
    function isotope_nms(a::Orbital, b::Orbital, Z::Float64, grid::Radial.Grid)
        mtp = min(size(a.P, 1), size(b.P, 1));   lb = Basics.subshell_l(b.subshell);    jb2   = Basics.subshell_2j(b.subshell)
        alphaZ = Defaults.getDefaults("alpha") * Z

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            error("stop a")
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp
                wb = (a.Pprime[i] * b.Pprime[i]  +  a.Qprime[i] * b.Qprime[i])  +
                     (lb*(lb+1) * a.P[i] * b.P[i]  +  (jb2-1)*jb2 * a.Q[i] * b.Q[i]) / (grid.r[i]^2)
                wc = - 2 * alphaZ * (a.Q[i] * b.Pprime[i]  +  b.Q[i] * a.Pprime[i]) / grid.r[i]
                wd = - alphaZ * (b.subshell.kappa - 1) * (a.Q[i] * b.P[i]  +  b.Q[i] * a.P[i]) / (grid.r[i]^2)
                wa = wa + (wb + wc + wd) * grid.wr[i]
            end
            return( wa )
        else
            error("stop b")
        end
    end

    """
    `RadialIntegrals.isotope_smsB(a::Orbital, c::Orbital, Z::Float64, grid::Radial.Grid)`
        ... computes the specific mass shift radial integral int_o^infty ... A value::Float64 is returned.
    """
    function isotope_smsB(a::Orbital, c::Orbital, Z::Float64, grid::Radial.Grid)
        mtp = min(size(a.P, 1), size(c.P, 1));   alphaZ = Defaults.getDefaults("alpha") * Z
        minusa = Subshell(1, -a.subshell.kappa);    minusc = Subshell(1, -c.subshell.kappa)

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            error("stop a")
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp
                wb = (- a.Q[i] * c.P[i] * AngularMomentum.sigma_reduced_me(minusa, c.subshell)  +
                        c.Q[i] * a.P[i] * AngularMomentum.sigma_reduced_me(a.subshell, minusc)  )
                wa = wa - alphaZ / grid.r[i] * wb * grid.wr[i]
            end
            return( wa )
        else
            error("stop b")
        end
    end

    """
    `RadialIntegrals.isotope_smsC(a::Orbital, c::Orbital, Z::Float64, grid::Radial.Grid)`
        ... computes the specific mass shift radial integral int_o^infty ... A value::Float64 is returned.
    """
    function isotope_smsC(a::Orbital, c::Orbital, Z::Float64, grid::Radial.Grid)
        mtp = min(size(a.P, 1), size(c.P, 1));   alphaZ = Defaults.getDefaults("alpha") * Z

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            error("stop a")
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp
                wa = wa - alphaZ / grid.r[i] * (a.Q[i] * c.P[i] - c.Q[i] * a.P[i]) * grid.wr[i]
            end
            return( wa )
        else
            error("stop b")
        end
    end

    """
    `RadialIntegrals.nondiagonalD(pm::Int64, kappa::Int64, bspline1::Bsplines.Bspline, bspline2::Bsplines.Bspline, grid::Radial.Grid)`
        ... computes the (radial and non-diagonal) D_kappa^+/- integral two the bsplines, all defined on grid
            <bspline1| +/- d/dr + kappa/r | bspline2>. -- pm = +1/-1 provides the phase for taking the derivative.
    """
    function nondiagonalD(pm::Int64, kappa::Int64, bspline1::Bsplines.Bspline, bspline2::Bsplines.Bspline, grid::Radial.Grid)
        if  bspline1.upper <= bspline2.lower  ||  bspline2.upper <= bspline1.lower    return( 0. )   end
        mtp = min( bspline1.upper, bspline2.upper)
        n0  = max( bspline1.lower, bspline2.lower)

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp

            function f1(i::Int64)
                wa = pm * bspline1.bs[i] * bspline2.bp[i]
                return( wa )
            end
            function f2(i :: Int64)
                if  i == 1  wa = bspline1.bs[i] * kappa * bspline2.bs[i] / (0.3 * grid.r[2])
                else        wa = bspline1.bs[i] * kappa * bspline2.bs[i] / grid.r[i]   end
                return( wa )
            end

            I1 = Math.integrateTransform(f1, n0, mtp, grid)
            I2 = Math.integrateTransform(f2, n0, mtp, grid)
            return( I1+I2 )
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = n0:mtp
                wa = wa + pm * bspline1.bs[i] * bspline2.bp[i] * grid.wr[i]
                if  i == 1  wa = wa + bspline1.bs[i] * kappa * bspline2.bs[i] / (0.3 * grid.r[2]) * grid.wr[i]
                else        wa = wa + bspline1.bs[i] * kappa * bspline2.bs[i] / grid.r[i] * grid.wr[i]   end
            end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.overlap()`

    + (orbital1::Radial.Orbital, orbital2::Radial.Orbital, grid::Radial.Grid)`
        ... computes the (radial) overlap integral <orbital_a|orbital_b>  for two relativistic orbitals of the same
            symmetry (kappa).
    """
    function overlap(orbital1::Radial.Orbital, orbital2::Radial.Orbital, grid::Radial.Grid)
        mtp = min(size(orbital1.P, 1), size(orbital2.P, 1))

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp

            function f(i :: Int64)
                return( orbital1.P[i] * orbital2.P[i] + orbital1.Q[i] * orbital2.Q[i] )
            end

            return( Math.integrateFitTransform(f, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 1:grid.nr
                if i > mtp   break   end
                wa = wa + ( orbital1.P[i] * orbital2.P[i] + orbital1.Q[i] * orbital2.Q[i] ) * grid.wr[i]
            end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    + (bspline1::Bsplines.Bspline, bspline2::Bsplines.Bspline, grid::Radial.Grid)`
        ... computes the (radial) overlap integral <bspline1|bsplines>  for two bpslines as defined on grid.
    """
    function overlap(bspline1::Bsplines.Bspline, bspline2::Bsplines.Bspline, grid::Radial.Grid)

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == Basics.MeshGrasp
            if  bspline1.upper <= bspline2.lower  ||  bspline2.upper <= bspline1.lower    return( 0. )   end
            mtp = min( bspline1.upper, bspline2.upper)
            n0  = max( bspline1.lower, bspline2.lower)

            function f(i :: Int64)
                return( bspline1.bs[i] * bspline2.bs[i] )
            end

            return( Math.integrateTransform(f, n0, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 1:grid.nr   wa = wa + bspline1.bs[i] * bspline2.bs[i] * grid.wr[i]   end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    + (p1List::Array{Float64,1}, p2List::Array{Float64,1}, grid::Radial.Grid)`
        ... computes the (radial) overlap integral of two (non-relativistic) radial orbital functions <p1|p2>  as defined on grid.
    """
    function overlap(p1List::Array{Float64,1}, p2List::Array{Float64,1}, grid::Radial.Grid)

        mtp = min( length(p1List), length(p2List))

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp

            function f(i :: Int64)
                return( p1List[i] * p2List[i] )
            end

            return( Math.integrateTransform(f, 1, mtp, grid) )
        elseif  grid.mesh == MeshGL
            error("stop a")
        else
            error("stop b")
        end
    end

    """
    `RadialIntegrals.qedDampedOverlap(lambda::Float64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)`
        ... computes the damped (radial) integral  int_0^infty (P_a P_b  +  Q_a Q_b) * e^{r/lambda} for the radial
            orbitals a, b on the given grid. A value::Float64 is returned.
    """
    function qedDampedOverlap(lambda::Float64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)
        mtp = min(size(a.P, 1), size(b.P, 1))
        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp   wb = Base.MathConstants.e^(- grid.r[i]/lambda);     wa = wa + (a.P[i]*wb*b.P[i] + a.Q[i]*wb*b.Q[i]) * grid.wr[i]   end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.qedLowFrequency(a::Radial.Orbital, b::Radial.Orbital, nm::Nuclear.Model, grid::Radial.Grid, qgrid::Radial.GridGL)`
        ... computes the (radial) integral for the low-frequency QED potential for the radial orbitals a, b on the given grid.
            A value::Float64 is returned.
    """
    function qedLowFrequency(a::Radial.Orbital, b::Radial.Orbital, nm::Nuclear.Model, grid::Radial.Grid, qgrid::Radial.GridGL)
        alpha = Defaults.getDefaults("alpha");    BZ = 0.074 + 0.035 * nm.Z * alpha
        mtp = min(size(a.P, 1), size(b.P, 1))
        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGL
            wa = 0.
            for  i = 1:mtp   wb = Base.MathConstants.e^(-nm.Z * grid.r[i]) ;     wa = wa + (a.P[i]*wb*b.P[i] + a.Q[i]*wb*b.Q[i]) * grid.wr[i]   end
            wa = -BZ * nm.Z^4 * alpha^3 * wa
         else
            error("stop a")
        end

        println("QED single-electron strength <$(a.subshell)| h^(SE, low-frequency) | $(b.subshell)> = $wa ")
        return( wa )
    end

    """
    `RadialIntegrals.qedUehling(a::Radial.Orbital, b::Radial.Orbital, pot::Radial.Potential,
                                grid::Radial.Grid, qgrid::Radial.GridGL)`
        ... computes the (radial) integral for the Uehling potential for the radial orbitals a, b on the given grid. This included a
            formal t-integration that is performed internally on the (QED) grid qgrid. A value::Float64 is returned.
    """
    function qedUehling(a::Radial.Orbital, b::Radial.Orbital, pot::Radial.Potential, grid::Radial.Grid, qgrid::Radial.GridGL)
        # Define the internal t-integration that is specific to the (simplified) Uehling potential; cf. PRA 72, 052115 (2005); eq. (9)
        function tIntegral(r::Float64)
            wx = 0.;
            for  i = 1:qgrid.nt   t = qgrid.t[i];    wx = wx + sqrt(t^2 - 1) / t^2 * (1 + 1 / (2*t^2)) * Base.MathConstants.e^(-2t*r) * qgrid.wt[i]   end
            return( wx )
        end

        mtp = min(size(a.P, 1), size(b.P, 1))
        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGL
            wa = 0.
            for  i = 1:mtp
                wb = tIntegral(grid.r[i]) * (-pot.Zr[i] * grid.r[i])
                wa = wa + (a.P[i]*wb*b.P[i] + a.Q[i]*wb*b.Q[i]) * grid.wr[i]
            end
            wa = 2. * Defaults.getDefaults("alpha") / (3pi) * wa
        else
            error("stop a")
        end

        println("QED single-electron strength <$(a.subshell)| h^(Uehling) | $(b.subshell)> = $wa ")
        return( wa )
    end

    """
    `RadialIntegrals.qedUehlingSimple(a::Radial.Orbital, b::Radial.Orbital, pot::Radial.Potential,
                                      grid::Radial.Grid, qgrid::Radial.GridGL)`
        ... computes the (radial) integral for the Uehling potential for the radial orbitals a, b on the given grid. This
            included a formal t-integration that is performed internally on the (QED) grid qgrid. A value::Float64 is returned.
    """
    function qedUehlingSimple(a::Radial.Orbital, b::Radial.Orbital, pot::Radial.Potential, grid::Radial.Grid, qgrid::Radial.GridGL)
        # Define the internal t-integration that is specific to the (simplified) Uehling potential; cf. PRA 72, 052115 (2005); eq. (9)
        function tIntegral(r::Float64)
            wx = 0.;
            for  i = 1:qgrid.nt   t = qgrid.t[i];    wx = wx + sqrt(t^2 - 1) / t^2 * (1 + 1 / (2*t^2)) * Base.MathConstants.e^(-2t*r) * qgrid.wt[i]   end
            return( wx )
        end

        mtp = min(size(a.P, 1), size(b.P, 1))
        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGL
            wa = 0.
            for  i = 1:mtp
                wb = tIntegral(grid.r[i]) * (-pot.Zr[i] * grid.r[i])
                wa = wa + (a.P[i]*wb*b.P[i] + a.Q[i]*wb*b.Q[i]) * grid.wr[i]
            end
            wa = 2. * Defaults.getDefaults("alpha") / (3pi) * wa
        else
            error("stop a")
        end

        println("QED single-electron strength <$(a.subshell)| h^(simplified Uehling) | $(b.subshell)> = $wa ")
        return( wa )
    end

    """
    `RadialIntegrals.rkDiagonal()`   ... computes the (radial and diagonal) integral of r^k for two radial orbital functions.

    + (k::Int64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)`
        ... computes this integral for two relativistic orbitals:   < r^k >_ab = int_0^\\infty  dr  [P_a P_b + Q_a Q_b]  r^k
    """
    function rkDiagonal(k::Int64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)
        mtp = min(size(a.P, 1), size(b.P, 1))

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            function f(i :: Int64)    return( (a.P[i] * b.P[i] + a.Q[i] * b.Q[i]) * (grid.r[i]^k) )    end
            return( Math.integrateFitTransform(f, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wa = 0.
            if  k > -3   m0 = 2   else   m0 = 6   end    # Don't allow too small r-values
            for  i = m0:mtp   wa = wa + (a.P[i] * b.P[i] + a.Q[i] * b.Q[i]) * (grid.r[i]^k) * grid.wr[i]   end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    + (k::Int64, p1List::Array{Float64,1}, p2List::Array{Float64,1}, grid::Radial.Grid)`
        ... computes this integral for two non-relativistic orbitals:   < r^k >_ab = int_0^\\infty  dr  [P_a P_b]  r^k
    """
    function rkDiagonal(k::Int64, p1List::Array{Float64,1}, p2List::Array{Float64,1}, grid::Radial.Grid)
        mtp = min( length(p1List), length(p2List))

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp

            function f(i :: Int64)
                return( p1List[i] * p2List[i] * (grid.r[i]^k) )
            end

            return( Math.integrateTransform(f, 2, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp   wa = wa + p1List[i] * p2List[i] * (grid.r[i]^k) * grid.wr[i]   end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.rkNonDiagonal(k::Int64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)`
        ... computes the (radial and non-diagonal) integral of r^k for two relativistic orbitals:
            [ r^k ]_ab = int_0^\\infty dr [P_a Q_b + Q_a P_b] r^k
    """
    function rkNonDiagonal(k::Int64, a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)
        mtp = min(size(a.P, 1), size(b.P, 1))

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            function f(i :: Int64)    return( (a.P[i] * b.Q[i] + a.Q[i] * b.P[i]) * (grid.r[i]^k) )    end
            return( Math.integrateFitTransform(f, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp   wa = wa + (a.P[i] * b.Q[i] + a.Q[i] * b.P[i]) * (grid.r[i]^k) * grid.wr[i]   end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.SlaterRk_2dim(k::Int64, a::Radial.Orbital, b::Radial.Orbital, c::Radial.Orbital, d::Orbital, grid::Radial.Grid)`
        ... computes the (relativistic) Slater integral

         R^k (abcd) = int_0^infty dr int_0^infty ds (P_a P_c + Q_a Q_c) r_<^k / r_>^(k+1) (P_b P_d + Q_b Q_d)

         of rank k for the four orbitals a, b, c, d, and over the given grid by using an explicit 2-dimensional integration scheme; a
         value::Float64 is returned.
    """
    function SlaterRk_2dim(k::Int64, a::Radial.Orbital, b::Radial.Orbital, c::Radial.Orbital, d::Radial.Orbital, grid::Radial.Grid)
        function ul(r :: Float64, s :: Float64) :: Float64
            if     r <= s    return( r^k/s^(k+1) )
            elseif r > s     return( s^k/r^(k+1) )
            end
        end

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            function fs(r :: Int64, s :: Int64) :: Float64
                return( ul(grid.r[r], grid.r[s]) * ( b.P[s] * d.P[s] + b.Q[s] * d.Q[s] ) )
            end

            function f(r :: Int64) :: Float64
            function ff(i :: Int64) :: Float64    return( fs(r, i) )    end
            return( (a.P[r] * c.P[r] + a.Q[r] * c.Q[r] ) * Math.integrateFitTransform(ff, min(size(b.P, 1), size(d.P, 1)), grid) )
            end

            return Math.integrateFitTransform(f, min(size(a.P, 1), size(c.P, 1)), grid)
        elseif  grid.mesh == MeshGL
            mtp_ac = min(size(a.P, 1), size(c.P, 1));    mtp_bd = min(size(b.P, 1), size(d.P, 1))
            wa = 0.
            for  r = 2:mtp_ac
                for  s = 2:mtp_bd   wa = wa + (a.P[r] * c.P[r] + a.Q[r] * c.Q[r]) * ul(grid.r[r], grid.r[s]) *
                                              (b.P[s] * d.P[s] + b.Q[s] * d.Q[s]) * grid.wr[r] * grid.wr[s]   end
            end
            ## println("Test: SlaterRk_2dim(); wa = $wa")
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.SlaterRk_new(k::Int64, a::Radial.Orbital, b::Radial.Orbital, c::Radial.Orbital, d::Orbital, grid::Radial.Grid)`
        ... computes the (relativistic) Slater integral

         R^k (abcd) = int_0^infty dr int_0^infty ds (P_a P_c + Q_a Q_c) r_<^k / r_>^(k+1) (P_b P_d + Q_b Q_d)

         of rank k for the four orbitals a, b, c, d, and over the given grid by using a factorized integration scheme; a
         value::Float64 is returned.
    """
    function SlaterRk_new(k::Int64, a::Radial.Orbital, b::Radial.Orbital, c::Radial.Orbital, d::Radial.Orbital, grid::Radial.Grid)
        function ul(r :: Float64, s :: Float64) :: Float64
            if     r <= s    return( r^k/s^(k+1) )
            elseif r > s     return( s^k/r^(k+1) )
            end
        end

        function rLowUp()
            rtuple = Tuple{Int64,Int64}[]
            for  n = 2:grid.orderGL:100000000
               if   n > grid.nr   break   end
               push!( rtuple, (n, n+grid.orderGL-1))
            end
            return( rtuple )
        end

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            error("stop a")
        elseif  grid.mesh == MeshGL
            mtp_ac = min(size(a.P, 1), size(c.P, 1));    mtp_bd = min(size(b.P, 1), size(d.P, 1));  rtuple = rLowUp()
            wa = 0.
            for  (rlow, rup)  in  rtuple
                for  (slow, sup)  in  rtuple
                    war = 0.;    was = 0.
                    if  sup < rlow
                        for  r = rlow:rup   if r > mtp_ac   break   end
                            war = war + (a.P[r] * c.P[r] + a.Q[r] * c.Q[r]) / grid.r[r]^(k+1) * grid.wr[r]    end
                        for  s = slow:sup   if s > mtp_bd   break   end
                            was = was + (b.P[s] * d.P[s] + b.Q[s] * d.Q[s]) * grid.r[s]^k     * grid.wr[s]    end
                        wa = wa + war * was
                    elseif  rup < slow
                        for  r = rlow:rup   if r > mtp_ac   break   end
                            war = war + (a.P[r] * c.P[r] + a.Q[r] * c.Q[r]) * grid.r[r]^k     * grid.wr[r]    end
                        for  s = slow:sup   if s > mtp_bd   break   end
                            was = was + (b.P[s] * d.P[s] + b.Q[s] * d.Q[s]) / grid.r[s]^(k+1) * grid.wr[s]    end
                        wa = wa + war * was
                    elseif  rlow == slow  &&  rup == sup
                        for  r = rlow:rup       if r > mtp_ac   break   end
                            for  s = slow:sup   if s > mtp_bd   break   end
                                wa = wa + (a.P[r] * c.P[r] + a.Q[r] * c.Q[r]) * ul(grid.r[r], grid.r[s]) *
                                          (b.P[s] * d.P[s] + b.Q[s] * d.Q[s]) * grid.wr[r] * grid.wr[s]       end
                        end
                    else  error("stop b")
                    end
                end
            end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.SlaterRk_DebyeHueckel_2dim(k::Int64, a::Radial.Orbital, b::Radial.Orbital, c::Radial.Orbital, d::Orbital,
                                                grid::Radial.Grid, lambda::Float64)`
        ... computes the (relativistic) Slater-Debye-Hueckel integral

            R^k (abcd) = int_0^infty dr int_0^infty ds (P_a P_c + Q_a Q_c) [r_<^k / r_>^(k+1)]^(DH screened) (P_b P_d + Q_b Q_d)

            of rank k for the four orbitals a, b, c, d, and over the given grid by using an explicit 2-dimensional integration
            scheme; a value::Float64 is returned.
    """
    function SlaterRk_DebyeHueckel_2dim(k::Int64, a::Radial.Orbital, b::Radial.Orbital, c::Radial.Orbital, d::Radial.Orbital,
                                        grid::Radial.Grid, lambda::Float64)

        function ul_DH(L::Int64, s::Float64, r::Float64)
            # Calculates the ul_DH(r,s) function for  s <= r.
            sum = 0.;   suma = 0.
            for  p = 0:2
                for q = 0:L
                    sum = sum + (2^(L-q)) * (lambda^(L+p+p-q)) *  factorial(L+q) * factorial(L+p) /
                                ( factorial(L+L+p+p+1) * factorial(L-q) * factorial(p) * factorial(q)) *
                                (s^(L+p+p)) * exp(-lambda*r) / (r^(q+1))
                end
                if (p == 2) suma = sum   end
            end
            return( (L+L+1) * sum )
        end

        function ul(r :: Float64, s :: Float64) :: Float64
            if     r <= s    return( ul_DH(k, r, s) )
            elseif r > s     return( ul_DH(k, s, r) )
            end
        end

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            function fs(r :: Int64, s :: Int64) :: Float64
                return( ul(grid.r[r], grid.r[s]) * ( b.P[s] * d.P[s] + b.Q[s] * d.Q[s] ) )
            end

            function f(r :: Int64) :: Float64
            function ff(i :: Int64) :: Float64    return( fs(r, i) )    end
            return( (a.P[r] * c.P[r] + a.Q[r] * c.Q[r] ) * Math.integrateFitTransform(ff, min(size(b.P, 1), size(d.P, 1)), grid) )
            end

            return Math.integrateFitTransform(f, min(size(a.P, 1), size(c.P, 1)), grid)
        elseif  grid.mesh == MeshGL
            mtp_ac = min(size(a.P, 1), size(c.P, 1));    mtp_bd = min(size(b.P, 1), size(d.P, 1))
            wa = 0.
            for  r = 2:mtp_ac
                for  s = 2:mtp_bd   wa = wa + (a.P[r] * c.P[r] + a.Q[r] * c.Q[r]) * ul(grid.r[r], grid.r[s]) *
                                              (b.P[s] * d.P[s] + b.Q[s] * d.Q[s]) * grid.wr[r] * grid.wr[s]   end
            end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.Vinti(a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)`
        ... computes the (radial) Vinti integral for the two radia integrals a and b:

                                            [ d        kappa_a (kappa_a+1) - kappa_b (kappa_b+1) ]
            R^(Vinti) = int_0^infty dr  P_a [ --   -  ------------------------------------------ ] P_b    +  similar (not equal) in Q_a, Q_b
                                            [ dr                          2 r                    ]

            a value::Float64 is returned.
    """
    function  Vinti(a::Radial.Orbital, b::Radial.Orbital, grid::Radial.Grid)

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            error("stop a")
        elseif  grid.mesh == MeshGL
            mtp_ab = min(size(a.P, 1), size(b.P, 1));    kapa = a.subshell.kappa;     kapb = b.subshell.kappa
            wa = 0.
            for  r = 2:mtp_ab
                wc = a.P[r] * b.Pprime[r] - a.P[r] * kapa * (kapa+1) * b.P[r] / (2grid.r[r])  + a.P[r] * kapb * (kapb+1) * b.P[r] / (2grid.r[r])
                wd = a.Q[r] * b.Qprime[r] + a.Q[r] * kapa * (-kapa+1)* b.Q[r] / (2grid.r[r]) - a.Q[r] * kapb * (-kapb+1) * b.Q[r] / (2grid.r[r])
                wa = wa +  (wc + wd) * grid.wr[r]
            end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.Vlocal(bspline1::Bsplines.Bspline, bspline2::Bsplines.Bspline, potential::Radial.Potential, grid::Radial.Grid)`
        ... computes the (radial) integral for the local potential and for two bpslines, all defined on grid
            <bspline1| potential.Zr | bspline2>. -- Here, potential. V must provide the effective charge zz(r) = - V * r.
    """
    function Vlocal(bspline1::Bsplines.Bspline, bspline2::Bsplines.Bspline, potential::Radial.Potential, grid::Radial.Grid)
        ## if  bspline1.upper <= bspline2.lower  ||  bspline2.upper <= bspline1.lower    return( 0. )   end
        mtp = min( bspline1.upper, bspline2.upper)
        n0  = max( bspline1.lower, bspline2.lower)

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp

            function f(i :: Int64)
                if  i == 1  wa = - bspline1.bs[i] * potential.Zr[i] * bspline2.bs[i] / (0.3 * grid.r[2])
                else        wa = - bspline1.bs[i] * potential.Zr[i] * bspline2.bs[i] / grid.r[i]   end
                return( wa )
            end

            return( Math.integrateTransform(f, n0, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = n0:mtp
                if  i == 1  wa = wa - bspline1.bs[i] * potential.Zr[i] * bspline2.bs[i] / (0.3 * grid.r[2]) * grid.wr[i]
                else        wa = wa - bspline1.bs[i] * potential.Zr[i] * bspline2.bs[i] / grid.r[i] * grid.wr[i]   end
            end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.V0(wa::Array{Float64,1}, mtp::Int64, grid::Radial.Grid)`
        ... computes the (radial) integral int_0^infty dr wa; a value::Float64 is returned.
    """
    function V0(wa::Array{Float64,1}, mtp::Int64, grid::Radial.Grid)

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            function f(i :: Int64)    return( wa[i] )     end
            return( Math.integrateFitTransform(f, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wb = 0.
            for  i = 1:mtp   wb = wb + wa[i] * grid.wr[i]   end
            return( wb )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.W5_Integral(mu::Int64, nu::Int64, orbitala::Radial.Orbital, orbitalc::Radial.Orbital,
                                                       orbitalb::Radial.Orbital, orbitald::Radial.Orbital, grid::Radial.Grid)`
        ... computes the (radial) integral for four relativistic orbitals:

            W_5 [ac|bd] = int_0^infty dr   int_0^r ds   [Pa Qc]_{r}  * ( s^nu / r^(nu+1) ) * [Pb Qd]_{s}

            as it frequently occurs in the frequency-independent Breit interaction.
    """
    function W5_Integral(mu::Int64, nu::Int64, a::Radial.Orbital, b::Radial.Orbital,
                                               c::Radial.Orbital, d::Radial.Orbital, grid::Radial.Grid)
        # Note mu = 5 is fixed historically and not used for this integral.
        !(mu == 5)   &&   error("mu = 5 required.")
        mtp = min(size(b.P, 1), size(d.P, 1))

        # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            function fs(s :: Int64) :: Float64
                return( grid.r[s]^nu * ( b.P[s] * d.Q[s] ) )
            end

            function f(r :: Int64) :: Float64
            function ff(i :: Int64) :: Float64  return( Math.integrateFitTransform(fs, i, grid) )  end
                if  r > mtp   return( 0. )
                else          return(  a.P[r] * c.Q[r] * ff(r) / grid.r[r]^(nu+1) )
                end
            end

            return Math.integrateFitTransform(f, min(size(a.P, 1), size(c.P, 1)), grid)
        elseif  grid.mesh == MeshGL
            mtp_ac = min(size(a.P, 1), size(c.P, 1));    mtp_bd = min(size(b.P, 1), size(d.P, 1))
            wa = 0.
            for  r = 2:mtp_ac
                for  s = 2:mtp_bd
                    if     s > r  continue
                    elseif s ==r  wa = wa + (a.P[r] * c.Q[r]) * (grid.r[s]^nu) / (grid.r[r]^(nu+1)) * (b.P[s] * d.Q[s]) * grid.wr[r] * grid.wr[s] / 2.0
                    else          wa = wa + (a.P[r] * c.Q[r]) * (grid.r[s]^nu) / (grid.r[r]^(nu+1)) * (b.P[s] * d.Q[s]) * grid.wr[r] * grid.wr[s]
                    end
                end
            end
            return( wa )
        else
            error("stop a")
        end
    end

    """
    `RadialIntegrals.Yk_ab(k::Int64, r::Float64, rho_ab::Array{Float64,1}, mtp::Int64, grid::Radial.Grid)`
        ... computes the (radial) integral

                                               r<^k
            Y_ab^k (r) = r * int_0^infty dr'  ------   rho_ab(r')
                                              r>^k+1

            a value::Float64 is returned.
    """
    function Yk_ab(k::Int64, r::Float64, rho_ab::Array{Float64,1}, mtp::Int64, grid::Radial.Grid)

       # Distinguish the radial integration for different grid definitions
        if  grid.mesh == MeshGrasp
            function f(i :: Int64)    rl = min(r, grid.r[i]);   rg = max(r, grid.r[i]);   return( rho_ab[i] * rl^k / rg^(k+1) )     end
            return( r * Math.integrateFitTransform(f, mtp, grid) )
        elseif  grid.mesh == MeshGL
            wa = 0.
            for  i = 2:mtp
                rl = min(r, grid.r[i]);   rg = max(r, grid.r[i])
                wa = wa + rho_ab[i] * rl^k / rg^(k+1) * grid.wr[i]
            end
            return( r * wa )
        else
            error("stop a")
        end
    end

end # module


text:
export DiracDelta, KroneckerDelta, DiracSpace, PointSpace

for TYP in (:DiracSpace,:PointSpace)
    @eval begin
        struct $TYP{T,D,R} <: Space{D,R}
          points::Vector{T}
          $TYP{T,D,R}(pts::AbstractVector{T}) where {T,D,R} = new(sort(pts))
        end

        function $TYP(points::AbstractVector{T}) where T
            $TYP{eltype(points),UnionDomain{Vector{Point{T}},Point{T}},real(T)}(points)
        end

        $TYP(points::Tuple) = $TYP(collect(points))
        $TYP() = $TYP(Float64[])
        $TYP(point::Number) = $TYP([point])
        $TYP(p::Point)=$TYP(p.x)

        dimension(d::$TYP)=length(d.points)

        # all points are equal, so only one block
        blocklengths(C::$TYP) = [length(C.points)]
        block(C::$TYP,k) = Block(1)

        domain(DS::$TYP) = UnionDomain(Point.(DS.points))
        setdomain(DS::$TYP,d::UnionDomain) = $TYP(map(d->d.x,components(d)))
        points(sp::$TYP,n::Integer)=sp.points[1:n]

        spacescompatible(a::$TYP,b::$TYP) = a.points == b.points
        canonicalspace(a::$TYP) = a

        union_rule(a::$TYP,b::$TYP) = $TYP(sort(union(a.points,b.points)))

        function coefficients(cfs::AbstractVector,fromspace::$TYP,tospace::$TYP)
            if spacescompatible(fromspace,tospace)
                return cfs
            end

            @assert length(cfs) ≤ length(fromspace.points)

            # this first for-loop removes coefficients of Dirac points that are zero
            nonzerofromspacepoints = eltype(fromspace.points)[]
            nonzeroDiraccfs = eltype(cfs)[]
            for i = 1:length(cfs)
                if cfs[i] != 0
                    push!(nonzerofromspacepoints, fromspace.points[i])
                    push!(nonzeroDiraccfs, cfs[i])
                end
            end

            # if the points that remain can be represented in the tospace
            if issubset(nonzerofromspacepoints,tospace.points)
                finalDiraccfs = eltype(cfs)[]
                j=1 #counter for the nonzerofromspacepoints
                for i = 1:length(tospace.points)
                    if j > length(nonzerofromspacepoints)
                        break
                    elseif nonzerofromspacepoints[j]==tospace.points[i]
                        push!(finalDiraccfs,nonzeroDiraccfs[j])
                        j += 1
                    else
                        push!(finalDiraccfs,0)
                    end
                end
                finalDiraccfs
            else
                error("The domain of the space being converted from has points that cannot be represented in the space you are converting to.")
            end
        end
    end
end

Space(d::Point) = PointSpace(d)

Fun(::typeof(identity), S::PointSpace) = Fun(S,S.points)
Fun(::typeof(identity), S::DiracSpace) = Fun(PointSpace(S.points),S.points)
transform(S::PointSpace,v::AbstractVector,plan...) = v
values(f::Fun{S}) where S<:PointSpace = coefficient(f,:)

function evaluate(f::AbstractVector,PS::PointSpace,x::Number)
    p = findfirst(y->isapprox(x,y),PS.points)
    if p == 0 || p > length(f)
        zero(eltype(f))
    else
        f[p]
    end
end

function evaluate(f::AbstractVector, PS::DiracSpace, x::Number)
    x ∉ domain(PS) && return zero(eltype(f))

    p = findfirst(y->isapprox(x,y), PS.points)
    if p == 0 || p > length(f)
        zero(eltype(f))
    else
        f[p]*Inf
    end
end

Base.sum(f::Fun{DS}) where DS<:DiracSpace =
    sum(f.coefficients[1:dimension(space(f))])

DiracDelta(x::Number)=Fun(DiracSpace(x),[1.])
DiracDelta()=DiracDelta(0.)

KroneckerDelta(x::Number) = Fun(PointSpace(x),[1.])
KroneckerDelta() = KroneckerDelta(0.)

function Base.cumsum(f::Fun{S},d::IntervalOrSegment{T}) where {S<:DiracSpace,T<:Real}
    pts=space(f).points
    @assert pts ==sort(pts)
    cfs=cumsum(f.coefficients)
    if leftendpoint(d) < first(pts) && rightendpoint(d) > last(pts)
        Fun(HeavisideSpace([leftendpoint(d);pts;rightendpoint(d)]),[0;cfs])
    elseif leftendpoint(d) == first(pts) && rightendpoint(d) > last(pts)
        Fun(HeavisideSpace([pts;rightendpoint(d)]),cfs)
    elseif leftendpoint(d) < first(pts) && rightendpoint(d) == last(pts)
        Fun(HeavisideSpace([leftendpoint(d);pts]),[0;cfs])
    elseif leftendpoint(d) == first(pts) && rightendpoint(d) == last(pts)
        Fun(HeavisideSpace(pts),cfs[1:end-1])
    else
        error("Implement")
    end
end

# for TYP in (:ReSpace,:Space)
#   @eval begin
#     function coefficients(cfs::AbstractVector,fromspace::$TYP,tospace::DiracSpace)
#       [0*tospace.points;coefficients(cfs,fromspace,tospace.space)]
#     end
#   end
# end

# function coefficients(cfs::AbstractVector,fromspace::DiracSpace,tospace::Space)
#   n = length(fromspace.points)
#   if n == 0 || cfs[1:n] == 0*cfs[1:n]
#       coefficients(fromspace.space,tospace)
#   else
#     error("The space you are converting from has Dirac deltas that cannot be represented in the space you are converting to.")
#   end
# end

function Multiplication(f::Fun{PS}, PS2::PointSpace) where PS<:PointSpace
    @assert space(f).points==PS2.points
    FiniteOperator(Diagonal(values(f)),PS2,PS2)
end

function Multiplication(f::Fun{PS}, DS::DiracSpace) where PS<:PointSpace
    @assert space(f).points==DS.points
    FiniteOperator(Diagonal(values(f)),DS,DS)
end

function Multiplication(f::Fun{DS}, PS::PointSpace) where DS<:DiracSpace
    @assert space(f).points==PS.points
    FiniteOperator(Diagonal(coefficient(f,:)),PS,space(f))
end

function coefficienttimes(f::Fun{PS}, g::Fun{DS}) where {PS<:PointSpace,DS<:DiracSpace}
    @assert space(f).points==space(g).points
    Fun(space(g),f.coefficients.*g.coefficients)
end

function coefficienttimes(f::Fun{DS}, g::Fun{PS}) where {PS<:PointSpace,DS<:DiracSpace}
    @assert space(f).points==space(g).points
    Fun(space(f),f.coefficients.*g.coefficients)
end

function coefficienttimes(f::Fun{PS}, g::Fun{PS2}) where {PS<:PointSpace,PS2<:PointSpace}
    @assert space(f).points==space(g).points
    Fun(space(g),f.coefficients.*g.coefficients)
end

/(f::Fun,g::Fun{PS}) where PS<:PointSpace = f*inv(g)

inv(f::Fun{PS}) where PS<:PointSpace = Fun(space(f),1 ./ f.coefficients)

#DiracSpace sampling
function randweights(pts, cfs)
    cs = cumsum(cfs)
    r = rand()
    if r≤cs[1]
        return pts[1]
    else
        for n=1:length(cfs)-1
            if cs[n]<r≤cs[n+1] && return pts[n+1]
                break
            end
        end
    end
end

function sample(f::Fun{<:DiracSpace,<:Real})
    cfs = f.coefficients/sum(f.coefficients)
    any(c -> c < 0, cfs) && throw(ArgumentError("All weights must be non-negative"))
    randweights(f.space.points, cfs)
end

function sample(f::Fun{<:DiracSpace,<:Real},n::Integer)
    p=[]
    for i=1:n
        p=[p;sample(f)]
    end
    return p
end
#integrate DiracSpace
function integrate(f::Fun{<:DiracSpace})
    pts=f.space.points
    cfs=f.coefficients
    int=Fun(HeavisideSpace([pts;Inf]),cumsum(cfs))
    return int
end

ext: jl


text:
## types for CM
# abstract
abstract type CM_Entity end
abstract type CM_Cell <: CM_Entity end
abstract type CM_Actual_Sum <: CM_Entity end
abstract type CM_Prediction_Sum <: CM_Entity end
abstract type CM_All_Sum <: CM_Entity end

# concrete
struct CM_TP <: CM_Cell end
struct CM_FP <: CM_Cell end
struct CM_FN <: CM_Cell end
struct CM_TN <: CM_Cell end

struct CM_AP <: CM_Actual_Sum end
struct CM_AN <: CM_Actual_Sum end
struct CM_PP <: CM_Prediction_Sum end
struct CM_PN <: CM_Prediction_Sum end
struct CM_ALL <: CM_All_Sum end

## Confusion Matrix
struct ConfusionMatrix
    tp::CM_TP
    fp::CM_FP
    fn::CM_FN
    tn::CM_TN
    ap::CM_AP
    an::CM_AN
    pp::CM_PP
    pn::CM_PN
    all::CM_ALL

    ConfusionMatrix() = new(CM_TP(), CM_FP(), CM_FN(), CM_TN(),
        CM_AP(), CM_AN(), CM_PP(), CM_PN(), CM_ALL())
end

struct CM_Value{T<:Number}
    tp::T
    fp::T
    fn::T
    tn::T
    ap::T
    an::T
    pp::T
    pn::T
    all::T

    function CM_Value(yhat::AbstractVector, y::AbstractVector)
        all = length(y)
        tp = yhat' * y
        ap = sum(y .== 1)
        pp = sum(yhat .== 1)

        an = all - ap
        pn = all - pp

        fp = pp - tp
        fn = ap - tp
        tn = an - fp

        return new{eltype(yhat)}(tp, fp, fn, tn, ap, an, pp, pn, all)
    end
end

# for storing constants
mutable struct ConstantOverPQ
    cPQ::Union{Matrix{Float32},Nothing}
    cPQ0::Union{Matrix{Float32},Nothing}
    cP0Q::Union{Matrix{Float32},Nothing}
    cP0Q0::Union{Matrix{Float32},Nothing}
    c::Union{Matrix{Float32},Nothing}

    ConstantOverPQ() = new(
        nothing, nothing, nothing, nothing, nothing)
end

function is_constant(co::ConstantOverPQ)
    return (isnothing(co.cPQ) && isnothing(co.cPQ0) && isnothing(co.cP0Q) &&
        isnothing(co.cP0Q0) && (!isnothing(co.c)))
end

function all_nothing(co::ConstantOverPQ)
    return (isnothing(co.cPQ) && isnothing(co.cPQ0) && isnothing(co.cP0Q) &&
        isnothing(co.cP0Q0) && isnothing(co.c))
end

function apply(func::Function, co::ConstantOverPQ, other)
    result = ConstantOverPQ()
    result.cPQ = isnothing(co.cPQ) ? nothing : func(co.cPQ, other)
    result.cPQ0 = isnothing(co.cPQ0) ? nothing : func(co.cPQ0, other)
    result.cP0Q = isnothing(co.cP0Q) ? nothing : func(co.cP0Q, other)
    result.cP0Q0 = isnothing(co.cP0Q0) ? nothing : func(co.cP0Q0, other)
    result.c = isnothing(co.c) ? nothing : func(co.c, other)
    return result
end

element_mult(x, y) = x .* y
element_div(x, y) = x ./ y
Base.:*(co::ConstantOverPQ, other) = apply(element_mult, co, other)
Base.:*(other, co::ConstantOverPQ) = apply(element_mult, co, other)
Base.:/(co::ConstantOverPQ, other) = apply(element_div, co, other)

function add_nothing(x, y)
    if isnothing(x) && isnothing(y)
        return nothing
    elseif isnothing(x)
        return y
    elseif isnothing(y)
        return x
    else
        return x + y
    end
end

function Base.:+(co::ConstantOverPQ, other::ConstantOverPQ)
    result = ConstantOverPQ()
    result.cPQ = add_nothing(co.cPQ, other.cPQ)
    result.cPQ0 = add_nothing(co.cPQ0, other.cPQ0)
    result.cP0Q = add_nothing(co.cP0Q, other.cP0Q)
    result.cP0Q0 = add_nothing(co.cP0Q0, other.cP0Q0)
    result.c = add_nothing(co.c, other.c)
    return result
end

## types for expression
abstract type CM_Expression end

# type for expr info
mutable struct ExpressionInfo
    is_linear_tp_tn::Bool       # is it linear w.r.t TP and TN
    depends_cell_cm::Bool       # does it depend on the cell of confusion matrix: tp, tn, fp, fn
    depends_actual_sum::Bool    # does it depend on actual sum statistics: AP & AN
    depends_predicted_sum::Bool # does it depend on predicted sum statistics: PP & PN
    is_constant::Bool           # does it contain only numbers or {ALL} entity
    needs_adv_sum_marg::Bool    # does it need 'sum'-marginal of the adversary to compute
    needs_pred_sum_marg::Bool   # does it need 'sum'-marginal of the predictor to compute
    is_constraint::Bool         # is it a constraint (contains '>=')

    # default values
    ExpressionInfo() = new(true, false, false, false, false, false, false, false)
end

# type for metric info
mutable struct MetricInfo
    valid::Bool                 # is it a valid metric based on our definition
    needs_adv_sum_marg::Bool    # does it need 'sum'-marginal of the adversary to compute
    needs_pred_sum_marg::Bool   # does it need 'sum'-marginal of the predictor to compute
    special_case_positive::Bool    # does it enforce special case for yhat = \zerovec && y = \zerovec
    special_case_negative::Bool   # does it enforce special case for yhat = \onevec && y = \onevec
    n_constraints::Int         # how many constraints does it have

    cs_special_case_positive_list::Vector{Bool}     # special cases for constraints
    cs_special_case_negative_list::Vector{Bool}     # special cases for constraints

    # default values
    MetricInfo() = new(false, false, false, false, false, 0, Bool[], Bool[])
end

# concrete expr
struct EXPR_UnaryEntity{TE <: CM_Entity} <: CM_Expression
    entity::TE
    multiplier::Float32
    info::ExpressionInfo

    function EXPR_UnaryEntity(entity::CM_Cell, multiplier = 1f0)
        info = ExpressionInfo()
        info.depends_cell_cm = true
        return new{typeof(entity)}(entity, multiplier, info)
    end

    function EXPR_UnaryEntity(entity::CM_Actual_Sum, multiplier = 1f0)
        info = ExpressionInfo()
        info.depends_actual_sum = true
        return new{typeof(entity)}(entity, multiplier, info)
    end

    function EXPR_UnaryEntity(entity::CM_Prediction_Sum, multiplier = 1f0)
        info = ExpressionInfo()
        info.depends_predicted_sum = true
        return new{typeof(entity)}(entity, multiplier, info)
    end

    function EXPR_UnaryEntity(entity::CM_All_Sum, multiplier = 1f0)
        info = ExpressionInfo()
        info.is_constant = true
        return new{typeof(entity)}(entity, multiplier, info)
    end
end

struct EXPR_UnaryIdentity <: CM_Expression
    multiplier::Float32
    info::ExpressionInfo

    function EXPR_UnaryIdentity(multiplier = 1f0)
        info = ExpressionInfo()
        info.is_constant = true
        return new(multiplier, info)
    end
end

struct EXPR_UnaryExpr{TE <: CM_Expression} <: CM_Expression
    expr::TE
    multiplier::Float32
    info::ExpressionInfo

    function EXPR_UnaryExpr(expr, multiplier = 1f0)
        return new{typeof(expr)}(expr, multiplier, expr.info)
    end
end

struct EXPR_Fraction{T1 <: CM_Expression, T2 <: CM_Expression} <: CM_Expression
    numerator::T1
    denominator::T2
    info::ExpressionInfo

    function EXPR_Fraction(numerator, denominator)
        info = ExpressionInfo()

        if denominator.info.is_constant
            info.is_linear_tp_tn = numerator.info.is_linear_tp_tn
            info.depends_cell_cm = numerator.info.depends_cell_cm
            info.depends_actual_sum = numerator.info.depends_actual_sum
            info.depends_predicted_sum = numerator.info.depends_predicted_sum
            info.is_constant = numerator.info.is_constant
        else
            if denominator.info.depends_cell_cm
                info.is_linear_tp_tn = false
            else
                info.is_linear_tp_tn = numerator.info.is_linear_tp_tn
            end
            info.depends_cell_cm = numerator.info.depends_cell_cm || denominator.info.depends_cell_cm
            info.depends_actual_sum = numerator.info.depends_actual_sum || denominator.info.depends_actual_sum
            info.depends_predicted_sum = numerator.info.depends_predicted_sum || denominator.info.depends_predicted_sum
            info.is_constant = false
        end

        info.needs_adv_sum_marg = denominator.info.depends_actual_sum ? true : numerator.info.needs_adv_sum_marg
        info.needs_pred_sum_marg = denominator.info.depends_predicted_sum ? true : numerator.info.needs_pred_sum_marg
        info.is_constraint = numerator.info.is_constraint || denominator.info.is_constraint

        return new{typeof(numerator), typeof(denominator)}(numerator, denominator, info)
    end
end

struct EXPR_Multiplication{T1 <: CM_Expression, T2 <: CM_Expression} <: CM_Expression
    expr1::T1
    expr2::T2
    info::ExpressionInfo

    function EXPR_Multiplication(expr1, expr2)
        info = ExpressionInfo()

        if !expr1.info.depends_cell_cm
            info.is_linear_tp_tn = expr2.info.is_linear_tp_tn
        elseif !expr2.info.depends_cell_cm
            info.is_linear_tp_tn = expr1.info.is_linear_tp_tn
        else
            info.is_linear_tp_tn = false
        end
        info.depends_cell_cm = expr1.info.depends_cell_cm || expr2.info.depends_cell_cm
        info.depends_actual_sum = expr1.info.depends_actual_sum || expr2.info.depends_actual_sum
        info.depends_predicted_sum = expr1.info.depends_predicted_sum || expr2.info.depends_predicted_sum
        info.is_constant = expr1.info.is_constant && expr2.info.is_constant

        if expr1.info.depends_actual_sum
            if expr2.info.is_constant
                info.needs_adv_sum_marg = expr1.info.needs_adv_sum_marg
            else
                info.needs_adv_sum_marg = true
            end
        elseif expr2.info.depends_actual_sum
            if expr1.info.is_constant
                info.needs_adv_sum_marg = expr2.info.needs_adv_sum_marg
            else
                info.needs_adv_sum_marg = true
            end
        else
            info.needs_adv_sum_marg = false
        end

        if expr1.info.depends_predicted_sum
            if expr2.info.is_constant
                info.needs_pred_sum_marg = expr1.info.needs_pred_sum_marg
            else
                info.needs_pred_sum_marg = true
            end
        elseif expr2.info.depends_predicted_sum
            if expr1.info.is_constant
                info.needs_pred_sum_marg = expr2.info.needs_pred_sum_marg
            else
                info.needs_pred_sum_marg = true
            end
        else
            info.needs_pred_sum_marg = false
        end

        info.is_constraint = expr1.info.is_constraint || expr2.info.is_constraint

        return new{typeof(expr1), typeof(expr2)}(expr1, expr2, info)
    end
end

struct EXPR_Power{TE <: CM_Expression} <: CM_Expression
    expr::TE
    power::Float32
    info::ExpressionInfo

    function EXPR_Power(expr, power)
        info = ExpressionInfo()

        if expr.info.depends_cell_cm
            info.is_linear_tp_tn = false
        else
            info.is_linear_tp_tn = expr.info.is_linear_tp_tn
        end
        info.depends_cell_cm = expr.info.depends_cell_cm
        info.depends_actual_sum = expr.info.depends_actual_sum
        info.depends_predicted_sum = expr.info.depends_predicted_sum
        info.is_constant = expr.info.is_constant
        info.needs_adv_sum_marg = expr.info.depends_actual_sum ? true : expr.info.needs_adv_sum_marg
        info.needs_pred_sum_marg = expr.info.depends_predicted_sum ? true : expr.info.needs_pred_sum_marg
        info.is_constraint = expr.info.is_constraint

        return new{typeof(expr)}(expr, power, info)
    end
end

struct EXPR_FunctionCall{TE <: CM_Expression} <: CM_Expression
    func::Function
    expr::TE
    info::ExpressionInfo

    function EXPR_FunctionCall(func, expr)
        info = ExpressionInfo()

        if expr.info.depends_cell_cm
            info.is_linear_tp_tn = false
        else
            info.is_linear_tp_tn = expr.info.is_linear_tp_tn
        end
        info.depends_cell_cm = expr.info.depends_cell_cm
        info.depends_actual_sum = expr.info.depends_actual_sum
        info.depends_predicted_sum = expr.info.depends_predicted_sum
        info.is_constant = expr.info.is_constant
        info.needs_adv_sum_marg = expr.info.depends_actual_sum ? true : expr.info.needs_adv_sum_marg
        info.needs_pred_sum_marg = expr.info.depends_predicted_sum ? true : expr.info.needs_pred_sum_marg
        info.is_constraint = expr.info.is_constraint

        return new{typeof(expr)}(func, expr, info)
    end
end

struct EXPR_Addition{T1 <: CM_Expression, T2 <: CM_Expression} <: CM_Expression
    expr1::T1
    expr2::T2
    info::ExpressionInfo

    function EXPR_Addition(expr1, expr2)
        info = ExpressionInfo()
        info.is_linear_tp_tn = expr1.info.is_linear_tp_tn && expr2.info.is_linear_tp_tn
        info.depends_cell_cm = expr1.info.depends_cell_cm || expr2.info.depends_cell_cm
        info.depends_actual_sum = expr1.info.depends_actual_sum || expr2.info.depends_actual_sum
        info.depends_predicted_sum = expr1.info.depends_predicted_sum || expr2.info.depends_predicted_sum
        info.is_constant = expr1.info.is_constant && expr2.info.is_constant
        info.needs_adv_sum_marg = expr1.info.needs_adv_sum_marg || expr2.info.needs_adv_sum_marg
        info.needs_pred_sum_marg = expr1.info.needs_pred_sum_marg || expr2.info.needs_pred_sum_marg
        info.is_constraint = expr1.info.is_constraint || expr2.info.is_constraint

        return new{typeof(expr1), typeof(expr2)}(expr1, expr2, info)
    end
end

struct EXPR_Constraint{TE <: CM_Expression} <: CM_Expression
    expr::TE
    threshold::Float32
    info::ExpressionInfo

    function EXPR_Constraint(expr, threshold)
        info = ExpressionInfo()

        info.is_linear_tp_tn = false
        info.depends_cell_cm = expr.info.depends_cell_cm
        info.depends_actual_sum = expr.info.depends_actual_sum
        info.depends_predicted_sum = expr.info.depends_predicted_sum
        info.is_constant = expr.info.is_constant
        info.needs_adv_sum_marg = expr.info.depends_actual_sum ? true : expr.info.needs_adv_sum_marg
        info.needs_pred_sum_marg = expr.info.depends_predicted_sum ? true : expr.info.needs_pred_sum_marg
        info.is_constraint = true

        return new{typeof(expr)}(expr, threshold, info)
    end
end

### functions on entities & expressions

## Pos and neg
Base.:+(e::CM_Entity) = EXPR_UnaryEntity(e, 1f0)
Base.:-(e::CM_Entity) = EXPR_UnaryEntity(e, -1f0)

Base.:+(e::CM_Expression) = e
Base.:-(e::CM_Expression) = EXPR_UnaryExpr(e, -1f0)

Base.:-(e::EXPR_UnaryEntity) = EXPR_UnaryEntity(e.entity, -e.multiplier)
Base.:-(e::EXPR_UnaryIdentity) = EXPR_UnaryIdentity(-e.multiplier)
Base.:-(e::EXPR_UnaryExpr) = EXPR_UnaryExpr(e.expr, -e.multiplier)

## Addition
# special cases entity
Base.:+(e1::CM_TP, e2::CM_FP) = CM_PP()
Base.:+(e1::CM_FP, e2::CM_TP) = CM_PP()
Base.:+(e1::CM_TP, e2::CM_FN) = CM_AP()
Base.:+(e1::CM_FN, e2::CM_TP) = CM_AP()

Base.:+(e1::CM_TN, e2::CM_FN) = CM_PN()
Base.:+(e1::CM_FN, e2::CM_TN) = CM_PN()
Base.:+(e1::CM_TN, e2::CM_FP) = CM_AN()
Base.:+(e1::CM_FP, e2::CM_TN) = CM_AN()

Base.:+(e1::CM_AP, e2::CM_AN) = CM_ALL()
Base.:+(e1::CM_AN, e2::CM_AP) = CM_ALL()
Base.:+(e1::CM_PP, e2::CM_PN) = CM_ALL()
Base.:+(e1::CM_PN, e2::CM_PP) = CM_ALL()

# general
Base.:+(e1::CM_Expression, e2::CM_Expression) = EXPR_Addition(e1, e2)
Base.:+(e1::CM_Expression, e2::CM_Entity) = e1 + EXPR_UnaryEntity(e2)
Base.:+(e1::CM_Entity, e2::CM_Expression) = EXPR_UnaryEntity(e1) + e2
Base.:+(e1::CM_Entity, e2::CM_Entity) = EXPR_UnaryEntity(e1) + EXPR_UnaryEntity(e2)

Base.:+(e1::CM_Expression, x::Number) = e1 + EXPR_UnaryIdentity(x)
Base.:+(x::Number, e1::CM_Expression) = EXPR_UnaryIdentity(x) + e1
Base.:+(e::CM_Entity, x::Number) = EXPR_UnaryEntity(e) + EXPR_UnaryIdentity(x)
Base.:+(x::Number, e::CM_Entity) = EXPR_UnaryIdentity(x) + EXPR_UnaryEntity(e)

## Subtraction
Base.:-(e1::CM_Expression, e2::CM_Expression) = e1 + (-e2)
Base.:-(e1::CM_Expression, e2::CM_Entity) = e1 + (-EXPR_UnaryEntity(e2))
Base.:-(e1::CM_Entity, e2::CM_Expression) = EXPR_UnaryEntity(e1) + (-e2)
Base.:-(e1::CM_Entity, e2::CM_Entity) = EXPR_UnaryEntity(e1) + (-e2)

Base.:-(e1::CM_Expression, e2) = e1 + (-e2)
Base.:-(e1::CM_Entity, e2) = EXPR_UnaryEntity(e1) + (-e2)
Base.:-(e1, e2::CM_Expression) = e1 + (-e2)
Base.:-(e1, e2::CM_Entity) = e1 + (-EXPR_UnaryEntity(e2))

## Multiplication
Base.:*(e1::CM_Expression, e2::CM_Expression) = EXPR_Multiplication(e1, e2)
Base.:*(e1::CM_Expression, e2::CM_Entity) = EXPR_Multiplication(e1, EXPR_UnaryEntity(e2))
Base.:*(e1::CM_Entity, e2::CM_Expression) = EXPR_Multiplication(EXPR_UnaryEntity(e1), e2)
Base.:*(e1::CM_Entity, e2::CM_Entity) = EXPR_Multiplication(EXPR_UnaryEntity(e1), EXPR_UnaryEntity(e2))

Base.:*(e::CM_Expression, x::Number) = EXPR_UnaryExpr(e, x)
Base.:*(x::Number, e::CM_Expression) = e * x
Base.:*(e::CM_Entity, x::Number) = EXPR_UnaryEntity(e, x)
Base.:*(x::Number, e::CM_Entity) = e * x

Base.:*(e::EXPR_UnaryEntity, x::Number) = EXPR_UnaryExpr(e.entity, x * e.multiplier)
Base.:*(e::EXPR_UnaryIdentity, x::Number) = EXPR_UnaryIdentity(x * e.multiplier)
Base.:*(e::EXPR_UnaryExpr, x::Number) = EXPR_UnaryExpr(e.expr, x * e.multiplier)

## Division
Base.:/(e1::CM_Expression, e2::CM_Expression) = EXPR_Fraction(e1, e2)
Base.:/(e1::CM_Expression, e2::CM_Entity) = EXPR_Fraction(e1, EXPR_UnaryEntity(e2))
Base.:/(e1::CM_Entity, e2::CM_Expression) = EXPR_Fraction(EXPR_UnaryEntity(e1), e2)
Base.:/(e1::CM_Entity, e2::CM_Entity) = EXPR_Fraction(EXPR_UnaryEntity(e1), EXPR_UnaryEntity(e2))

Base.:/(e::CM_Expression, x::Number) = EXPR_UnaryExpr(e, 1/x)
Base.:/(x::Number, e::CM_Expression) = EXPR_UnaryIdentity(x) / e
Base.:/(e::CM_Entity, x::Number) = EXPR_UnaryEntity(e, 1/x)
Base.:/(x::Number, e::CM_Entity) = EXPR_UnaryIdentity(x) / e

Base.:/(e::EXPR_UnaryEntity, x::Number) = EXPR_UnaryExpr(e.entity, e.multiplier / x)
Base.:/(e::EXPR_UnaryExpr, x::Number) = EXPR_UnaryExpr(e.expr, e.multiplier / x)

Base.:/(e::EXPR_UnaryIdentity, x::Number) = EXPR_UnaryIdentity(e.multiplier / x)
Base.:/(x::Number, e::EXPR_UnaryIdentity) = EXPR_UnaryIdentity(x / e.multiplier)

## Power
Base.:^(e::CM_Expression, x::Number) = EXPR_Power(e, x)
Base.:^(e::CM_Entity, x::Number) = EXPR_Power(EXPR_UnaryEntity(e), x)

## Greater than equal (for metric constraints)
Base.:>=(e::CM_Expression, x::Number) = EXPR_Constraint(e, x)
Base.:>=(e::CM_Entity, x::Number) = EXPR_Constraint(EXPR_UnaryEntity(e), x)

## Function call
Base.sqrt(e::CM_Expression) = EXPR_FunctionCall(Base.sqrt, e)
Base.sqrt(e::CM_Entity) = EXPR_FunctionCall(Base.sqrt, EXPR_UnaryEntity(e))
Base.log(e::CM_Expression) = EXPR_FunctionCall(Base.log, e)
Base.log(e::CM_Entity) = EXPR_FunctionCall(Base.log, EXPR_UnaryEntity(e))
Base.exp(e::CM_Expression) = EXPR_FunctionCall(Base.exp, e)
Base.exp(e::CM_Entity) = EXPR_FunctionCall(Base.exp, EXPR_UnaryEntity(e))

## Printing
Base.show(io::IO, x::CM_TP) =  print(io, "TP")
Base.show(io::IO, x::CM_FP) =  print(io, "FP")
Base.show(io::IO, x::CM_FN) =  print(io, "FN")
Base.show(io::IO, x::CM_TN) =  print(io, "TN")

Base.show(io::IO, x::CM_AP) =  print(io, "AP")
Base.show(io::IO, x::CM_AN) =  print(io, "AN")
Base.show(io::IO, x::CM_PP) =  print(io, "PP")
Base.show(io::IO, x::CM_PN) =  print(io, "PN")
Base.show(io::IO, x::CM_ALL) =  print(io, "ALL")

Base.show(io::IO, x::CM_Entity) = print(io, SubString(string(typeof(x)), 4))
Base.show(io::IO, x::CM_Expression) = print(io, "Expression::" * typeof(x))
Base.show(io::IO, x::EXPR_UnaryIdentity) = print(io, x.multiplier)

function Base.show(io::IO, x::EXPR_UnaryEntity)
    if x.multiplier == 1
        print(io, x.entity)
    elseif x.multiplier == -1
        print(io, "-" * string(x.entity))
    else
        print(io, string(x.multiplier) * " " * string(x.entity))
    end
end

function Base.show(io::IO, x::EXPR_UnaryExpr)
    if x.multiplier == 1
        print(io, x.expr)
    elseif x.multiplier == -1
        print(io, "-" * "(" * string(x.expr) * ")")
    else
        print(io, string(x.multiplier) * " " * "(" * string(x.expr) * ")")
    end
end

Base.show(io::IO, x::EXPR_Fraction) = print(
    io, "(" * string(x.numerator) * ")" * " / " * "(" * string(x.denominator) * ")")
Base.show(io::IO, x::EXPR_Multiplication) = print(
        io, "(" * string(x.expr1) * ")" * " * " * "(" * string(x.expr2) * ")")
Base.show(io::IO, x::EXPR_Addition) = print(
    io, "(" * string(x.expr1) * ")" * " + " * "(" * string(x.expr2) * ")")

Base.show(io::IO, x::EXPR_Power) = print(io, "(" * string(x.expr) * ")^" * string(x.power))
Base.show(io::IO, x::EXPR_FunctionCall) = print(io, string(x.func) * "(" * string(x.expr) * ")" )
Base.show(io::IO, x::EXPR_Constraint) = print(io, "(" * string(x.expr) * " >= " * string(x.threshold))

## Compute value
compute_value(::CM_TP, Cval::CM_Value) = Cval.tp
compute_value(::CM_FP, Cval::CM_Value) = Cval.fp
compute_value(::CM_FN, Cval::CM_Value) = Cval.fn
compute_value(::CM_TN, Cval::CM_Value) = Cval.tn

compute_value(::CM_AP, Cval::CM_Value) = Cval.ap
compute_value(::CM_AN, Cval::CM_Value) = Cval.an
compute_value(::CM_PP, Cval::CM_Value) = Cval.pp
compute_value(::CM_PN, Cval::CM_Value) = Cval.pn
compute_value(::CM_ALL, Cval::CM_Value) = Cval.all

compute_value(e::EXPR_UnaryEntity, Cval::CM_Value) = e.multiplier * compute_value(e.entity, Cval)
compute_value(e::EXPR_UnaryIdentity, Cval::CM_Value) = e.multiplier
compute_value(e::EXPR_UnaryExpr, Cval::CM_Value) = e.multiplier * compute_value(e.expr, Cval)

compute_value(e::EXPR_Fraction, Cval::CM_Value) =
    compute_value(e.numerator, Cval) / compute_value(e.denominator, Cval)
compute_value(e::EXPR_Multiplication, Cval::CM_Value) =
    compute_value(e.expr1, Cval) * compute_value(e.expr2, Cval)
compute_value(e::EXPR_Addition, Cval::CM_Value) =
    compute_value(e.expr1, Cval) + compute_value(e.expr2, Cval)

compute_value(e::EXPR_Power, Cval::CM_Value) = compute_value(e.expr, Cval) ^ e.power
compute_value(e::EXPR_FunctionCall, Cval::CM_Value) = e.func(compute_value(e.expr, Cval))

compute_value(e::EXPR_Constraint, Cval::CM_Value) = compute_value(e.expr, Cval) >= e.threshold

compute_value(::Nothing, ::CM_Value) = 0f0

## Compute value, m = n + 1, for \sum y_i = 0 to n
function compute_scaling(::CM_TP, m::Integer, info::MetricInfo)
    res = ConstantOverPQ()
    res.cPQ = ones(Float32, m,m)
    return res
end
function compute_scaling(::CM_FP, m::Integer, info::MetricInfo)
    res = ConstantOverPQ()
    res.cPQ0 = ones(Float32, m,m)
    return res
end
function compute_scaling(::CM_FN, m::Integer, info::MetricInfo)
    res = ConstantOverPQ()
    res.cP0Q = ones(Float32, m,m)
    return res
end
function compute_scaling(::CM_TN, m::Integer, info::MetricInfo)
    res = ConstantOverPQ()
    res.cP0Q0 = ones(Float32, m,m)
    return res
end

function compute_scaling(::CM_AP, m::Integer, info::MetricInfo)
    res = ConstantOverPQ()
    n = m - 1
    res.c = Float32.(repeat(0:n, 1, m)')
    return res
end
function compute_scaling(::CM_AN, m::Integer, info::MetricInfo)
    res = ConstantOverPQ()
    n = m - 1
    res.c = Float32.(repeat(n:-1:0, 1, m)')
    return res
end
function compute_scaling(::CM_PP, m::Integer, info::MetricInfo)
    res = ConstantOverPQ()
    n = m - 1
    res.c = Float32.(repeat(0:n, 1, m))
    return res
end
function compute_scaling(::CM_PN, m::Integer, info::MetricInfo)
    res = ConstantOverPQ()
    n = m - 1
    res.c = Float32.(repeat(n:-1:0, 1, m))
    return res
end
function compute_scaling(::CM_ALL, m::Integer, info::MetricInfo)
    res = ConstantOverPQ()
    n = m - 1
    res.c = n * ones(Float32, m,m)
    return res
end

function compute_scaling(e::EXPR_UnaryEntity, m::Integer, info::MetricInfo)
    return e.multiplier * compute_scaling(e.entity, m, info)
end
function compute_scaling(e::EXPR_UnaryIdentity, m::Integer, info::MetricInfo)
    res = ConstantOverPQ()
    res.c = e.multiplier * ones(Float32, m,m)
    return res
end
function compute_scaling(e::EXPR_UnaryExpr, m::Integer, info::MetricInfo)
    return e.multiplier * compute_scaling(e.expr, m, info)
end

function compute_scaling(e::EXPR_Fraction, m::Integer, info::MetricInfo)
    res_num = compute_scaling(e.numerator, m, info)
    res_den = compute_scaling(e.denominator, m, info)

    if is_constant(res_den)
        C = res_den.c ::Matrix{Float32}
        invC = zeros(Float32, m,m)

        if info.special_case_positive && info.special_case_negative
            invC[2:m-1, 2:m-1] = 1 ./ C[2:m-1, 2:m-1]
        elseif info.special_case_positive
            invC[2:m, 2:m] = 1 ./ C[2:m, 2:m]
        elseif info.special_case_negative
            invC[1:m-1,1:m-1] = 1 ./ C[1:m-1,1:m-1]
        else
            invC = 1 ./ C
        end

        # handle diviion by zero, when special cases not set => set it to zero
        invC[isinf.(invC)] .= 0f0

        return res_num * invC
    else
        throw(ErrorException("The metric contains unsupported operations"))
    end
end
function compute_scaling(e::EXPR_Multiplication, m::Integer, info::MetricInfo)
    res_ex1 = compute_scaling(e.expr1, m, info)
    res_ex2 = compute_scaling(e.expr2, m, info)

    if is_constant(res_ex1)
        return res_ex2 * res_ex1.c
    elseif is_constant(res_ex2)
        return res_ex1 * res_ex2.c
    else
        throw(ErrorException("The metric contains unsupported operations"))
    end
end

function compute_scaling(e::EXPR_Addition, m::Integer, info::MetricInfo)
    return compute_scaling(e.expr1, m, info) + compute_scaling(e.expr2, m, info)
end

function compute_scaling(e::EXPR_Power, m::Integer, info::MetricInfo)
    res_expr = compute_scaling(e.expr, m, info)
    if is_constant(res_expr)
        res_expr.c = res_expr.c .^ e.power
        return res_expr
    else
        throw(ErrorException("The metric contains unsupported operations"))
    end
end
function compute_scaling(e::EXPR_FunctionCall, m::Integer, info::MetricInfo)
    res_expr = compute_scaling(e.expr, m, info)
    if is_constant(res_expr)
        res_expr.c = e.func.(res_expr.c)
        return res_expr
    else
        throw(ErrorException("The metric contains unsupported operations"))
    end
end

compute_scaling(::Nothing, ::Integer, ::MetricInfo) = ConstantOverPQ()


text:
using StatisticalRethinking

PrPV = 0.95
PrPM = 0.01
PrV = 0.001
PrP = PrPV*PrV + PrPM*(1-PrV)
PrVP = PrPV*PrV / PrP

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
abstract Auth

# -------

immutable APIKeyAuth <: Auth
    key
end

function show(io::IO, a::APIKeyAuth)
    key_str =  repeat("*", length(a.key) - 3) * a.key[(end - 2):end]
    print(io, "MapLight API Key ($key_str)")
end

function authentiate_query(query_params, a::APIKeyAuth)
    query_params["apikey"] = a.key
end

# -------

immutable AnonymousAuth <: Auth
end

function show(io::IO, a::AnonymousAuth)
    print(io, "Anonymous MapLight API Auth")
end

function authentiate_query(query_params, a::AnonymousAuth)
    query_params  # nothing to be done
end

# -------

function authenticate(key::String)
    APIKeyAuth(key)
end

function authenticate()
    AnonymousAuth()
end


text:
for i in 1:10
    println(i)
end
1+2

max_stars_repo_name: Aaron-Gregory/RE-QSVT-HHL

max_issues_repo_name: Aaron-Gregory/RE-QSVT-HHL

max_forks_repo_name: Aaron-Gregory/RE-QSVT-HHL


text:
#=
run(`$(ENV["HOME"])/local/bin/mcl`)
run(`~/local/mcl test.data --abc -o out.cathat -I 10`)
=#

using PyCall,DataFrames
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport ncdata
filename = "ropa_g_nhept.nc"

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
# # Serving

"""
Start a server to present built output of project `src` for the given
`targets`. The actions performed by each `target` are defined by the targets
themselves. The default target used is [`html`](#).
"""
function serve(src, targets...=html; port=8000, keywords...)
    @info "Serving '$src'..."
    old = pwd() # Save directory so it can be returned to in cleanup.
    p = Project(src; keywords...)
    dir = mktempdir()
    sw = LiveServer.SimpleWatcher()
    watch_files!(sw, p)
    callback = function (path) # Callback to update project and rerun
        @info "path: $path"    # given targets.
        revise(p)
        update!(p)
        @sync for f in targets
            @info "target: $f"
            @async f(p, dir)
        end
        watch_files!(sw, p)
    end
    LiveServer.set_callback!(sw, callback)
    LiveServer.start(sw)
    try
        ## Run target-specific initialisation code.
        @sync for f in targets
            @async init(p, f; port=port, dir=dir)
        end
    finally
        ## Cleanup steps.
        LiveServer.stop(sw)
        rm(dir; recursive=true, force=true)
        cd(old) # Reinstate original directory when server is terminated.
    end
end

"""
Set files to watch for changes within a given [`Project`](#) `p`.
"""
function watch_files!(sw::LiveServer.SimpleWatcher, p::Project)
    empty!(sw.watchedfiles) # TODO: proper API for this?
    for file in FileTrees.files(p.tree)
        f = string(path(file))
        LiveServer.is_watched(sw, f) || LiveServer.watch_file!(sw, f)
    end
end

"""
Function for defining the initialisation code called by [`serve`](#)
for a specific target output.

Should be defined for each valid output target, i.e [`html`](#), [`pdf`](#),
etc. Signature should follow the following template:

```julia
function init(project::Project, ::typeof(target); kws...)
    # ...
end
```

where `target` is the `Function` defining the output format.
"""
function init end


text:
import GLFW: Window, MouseButton, Action
"""
Standard set of callback functions
"""
function standard_callbacks()
    return Function[window_open, window_size, window_position, keyboard_buttons,
                    mouse_buttons, dropped_files, framebuffer_size, unicode_input,
                    cursor_position, scroll, hasfocus, entered_window,]
end

"""
Returns a signal, which is true as long as the window is open.
returns `Bool`
[GLFW Docs](http://www.glfw.org/docs/latest/group__window.html#gaade9264e79fae52bdb78e2df11ee8d6a)
"""
function window_open(window::Window, s::Observable{Bool} = Observable(false))
    GLFW.SetWindowCloseCallback(window, (window,) -> s[])
    return s
end

"""
Size of the window. Must not correlate to the real pixel size.
This is why there is also framebuffer_size.
returns `NTuple{2,Int}}`
[GLFW Docs](http://www.glfw.org/docs/latest/group__window.html#gaaca1c2715759d03da9834eac19323d4a)
"""
function window_size(window::Window,
                     s::Observable{NTuple{2,Int}} = Observable(Int.(values(GLFW.GetWindowSize(window)))))
    GLFW.SetWindowSizeCallback(window,
                               (window, w::Cint, h::Cint) -> begin
                                   s[] = (Int(w), Int(h))
                               end)
    return s
end
"""
Size of window in pixel.
returns `NTuple{2,Int}}`
[GLFW Docs](http://www.glfw.org/docs/latest/group__window.html#ga311bb32e578aa240b6464af494debffc)
"""
function framebuffer_size(window::Window,
                          s::Observable{NTuple{2,Int}} = Observable(Int.(values(GLFW.GetFramebufferSize(window)))))
    GLFW.SetFramebufferSizeCallback(window,
                                    (window, w::Cint, h::Cint) -> begin
                                        s[] = (Int(w), Int(h))
                                    end)
    return s
end
"""
Position of the window in screen coordinates.
returns `NTuple{2,Int}}`
[GLFW Docs](http://www.glfw.org/docs/latest/group__window.html#ga1c36e52549efd47790eb3f324da71924)
"""
function window_position(window::Window, s::Observable{NTuple{2,Int}} = Observable((0, 0)))
    GLFW.SetWindowPosCallback(window,
                              (window, x::Cint, y::Cint) -> begin
                                  s[] = (Int(x), Int(y))
                              end)
    return s
end
"""
Registers a callback for the mouse buttons + modifiers
returns `NTuple{4, Int}`
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#ga1e008c7a8751cea648c8f42cc91104cf)
"""
function keyboard_buttons(window::Window,
                          s::Observable{NTuple{4,Int}} = Observable((0, 0, 0, 0)))
    keydict = Dict{Int,Bool}()
    GLFW.SetKeyCallback(window,
                        (window, button::GLFW.Key, scancode::Cint, action::GLFW.Action, mods::Cint) -> begin
                            s[] = (Int(button), Int(scancode), Int(action), Int(mods))
                        end)
    return s
end
"""
Registers a callback for the mouse buttons + modifiers
returns an `NTuple{3, Int}`,
containing the pressed button the action and modifiers.
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#ga1e008c7a8751cea648c8f42cc91104cf)
"""
function mouse_buttons(window::Window, s::Observable{NTuple{3,Int}} = Observable((0, 0, 0)))
    GLFW.SetMouseButtonCallback(window,
                                (window, button::GLFW.MouseButton, action::GLFW.Action, mods::Cint) -> begin
                                    s[] = (Int(button), Int(action), Int(mods))
                                end)
    return s
end
"""
Registers a callback for drag and drop of files.
returns `Vector{String}`, which are absolute file paths
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#gacc95e259ad21d4f666faa6280d4018fd)
"""
function dropped_files(window::Window, s::Observable{Vector{String}} = Observable(String[]))
    GLFW.SetDropCallback(window, (window, files) -> begin
                             s = map(String, files)
                         end)
    return s
end

"""
Registers a callback for keyboard unicode input.
returns an `Vector{Char}`,
containing the pressed char. Is empty, if no key is pressed.
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#ga1e008c7a8751cea648c8f42cc91104cf)
"""
function unicode_input(window::Window, s::Observable{Vector{Char}} = Observable(Char[]))
    GLFW.SetCharCallback(window, (window, c::Char) -> begin
                             s = c
                         end)
    return s
end
"""
Registers a callback for the mouse cursor position.
returns an `NTuple{2, Float64}`,
which is not in screen coordinates, with the upper left window corner being 0
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#ga1e008c7a8751cea648c8f42cc91104cf)
"""
function cursor_position(window::Window,
                         s::Observable{NTuple{2,Float64}} = Observable((0.0, 0.0)))
    GLFW.SetCursorPosCallback(window, (window, x::Cdouble, y::Cdouble) -> begin
                                  s[] = (x, y)
                              end)
    return s
end
"""
Registers a callback for the mouse scroll.
returns an `NTuple{2, Float64}`,
which is an x and y offset.
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#gacc95e259ad21d4f666faa6280d4018fd)
"""
function scroll(window::Window, s::Observable{NTuple{2,Float64}} = Observable((0.0, 0.0)))
    GLFW.SetScrollCallback(window,
                           (window, xoffset::Cdouble, yoffset::Cdouble) -> begin
                               s[] = (s[][1] + xoffset, s[][2] + yoffset)
                           end)
    return s
end
"""
Registers a callback for the focus of a window.
returns a `Bool`,
which is true whenever the window has focus.
[GLFW Docs](http://www.glfw.org/docs/latest/group__window.html#ga6b5f973531ea91663ad707ba4f2ac104)
"""
function hasfocus(window::Window, s::Observable{Bool} = Observable(false))
    GLFW.SetWindowFocusCallback(window, (window, focus::Bool) -> begin
                                    s[] = focus
                                end)
    return s
end
"""
Registers a callback for if the mouse has entered the window.
returns a `Bool`,
which is true whenever the cursor enters the window.
[GLFW Docs](http://www.glfw.org/docs/latest/group__input.html#ga762d898d9b0241d7e3e3b767c6cf318f)
"""
function entered_window(window::Window, s::Observable{Bool} = Observable(false))
    GLFW.SetCursorEnterCallback(window, (window, entered::Bool) -> begin
                                    s[] = entered
                                end)
    return s
end

#Came from GLWindow/events.jl
"""
Builds a Set of keys, accounting for released and pressed keys
"""
function currently_pressed_keys(v0::Set{Int}, button_action_mods)
    button, action, mods = button_action_mods
    if button != GLFW.KEY_UNKNOWN
        if action == GLFW.PRESS
            push!(v0, button)
        elseif action == GLFW.RELEASE
            delete!(v0, button)
        elseif action == GLFW.REPEAT
            # nothing needs to be done, besides returning the same set of keys
        else
            @error "Unrecognized enum value for GLFW button press action: $action"
        end
    end
    return v0
end

function remove_scancode(button_scancode_action_mods)
    button, scancode, action, mods = button_scancode_action_mods
    return button, action, mods
end
isreleased(button) = button[2] == GLFW.RELEASE
isdown(button) = button[2] == GLFW.PRESS

#question: do we need this?
# """
# Creates high level signals from the raw GLFW button signals.
# Returns a dictionary with button released and down signals.
# It also creates a signal, which is the set of all currently pressed buttons.
# `name` is used to name the dictionary keys.
# `buttons` is a tuple of (button, action, mods)::NTuple{3, Int}
# """
# function button_signals(buttons::NTuple{3, Int}, name::Symbol)
#     keyset = Set{Int}()
#     sizehint!(keyset, 10) # make it less suspicable to growing/shrinking
#     released = filter(isreleased, buttons.value, buttons)
#     down     = filter(isdown, buttons.value, buttons)
#     Dict{Symbol, Any}(
#         Symbol("$(name)_released") => map(first, released),
#         Symbol("$(name)_down")     => map(first, down),
#         Symbol("$(name)s_pressed") => foldp(
#             currently_pressed_keys, keyset, buttons
#         )
#     )
# end


text:
"""
compute the kronecker product of a and b.\n
used as a shorthand notation.\n
"""
⊗(a, b) = kron(a, b)

"""
AToE = _electric_vector_potential_to_electric_field_operator(number_of_cells, resolution, k_b)\n
This operators is used to compute the total electric field on the FDTD grid (Yee cell) given the vector potential.\n
k_b is the wavenumber in the background medium.\n
"""
function _electric_vector_potential_to_electric_field_operator(
    number_of_cells,
    resolution,
    k_b
)
    T = ComplexF32
    p = number_of_cells[1]
    q = number_of_cells[2]
    r = number_of_cells[3]
    Δx = Float32(resolution[1])
    Δy = Float32(resolution[2])
    Δz = Float32(resolution[3])

    k_b = Float32(k_b)
    # needed -> Ex : (P+2) x (Q+1) : (R+1), wanted -> Ex : P x (Q+1) : (R+1)
    # needed -> Ey : (P+1) x (Q+2) : (R+1), wanted -> Ey : (P+1) x Q : (R+1)
    # needed -> Ez : (P+1) x (Q+1) : (R+2), wanted -> Ez : (P+1) x (Q+1) : R
    #operators for Ex
    Rp = spdiagm(p, p+2, 1=>ones(T, p))
    Rq = spdiagm(q+1, q+1, 0=>ones(T, q+1))
    Rr = spdiagm(r+1, r+1, 0=>ones(T, r+1))

    Ix   = (Rr ⊗ (Rq ⊗ Rp))
    ∂x   = 1/Δx .*spdiagm(p, p+1, 0 =>-one(T) .*ones(T, p), 1 =>ones(T, p))
    ∂y   = 1/Δy .*spdiagm(q+1, q+2, 0 =>-one(T) .*ones(T, q+1), 1 =>ones(T, q+1))
    ∂z   = 1/Δz .*spdiagm(r+1, r+2, 0 =>-one(T) .*ones(T, r+1), 1 =>ones(T, r+1))
    ∂x²  = 1/(Δx)^2 .*spdiagm(p, p+2, 0 => ones(T, p),1=> -2 .*ones(T, p), 2 => ones(T, p))

    xx = Ix'*((Rr ⊗ (Rq ⊗ ∂x²)) + k_b^2 .*Ix)
    yx = Ix'*((Rr ⊗ (∂y ⊗ ∂x)))
    zx = Ix'*((∂z ⊗ (Rq ⊗ ∂x)))

    #operators for Ey
    Rp = spdiagm(p+1, p+1, 0=>ones(T, p+1))
    Rq = spdiagm(q, q+2, 1=>ones(T, q))
    Rr = spdiagm(r+1, r+1, 0=>ones(T, r+1))

    Iy   = (Rr ⊗ (Rq ⊗ Rp))
    ∂x   = 1/Δx .*spdiagm(p+1, p+2, 0 =>-one(T) .*ones(T, p+1), 1 =>ones(T, p+1))
    ∂y   = 1/Δy .*spdiagm(q, q+1, 0 =>-one(T) .*ones(T, q), 1 =>ones(T, q))
    ∂z   = 1/Δz .*spdiagm(r+1, r+2, 0 =>-one(T) .*ones(T, r+1), 1 =>ones(T, r+1))
    ∂y²  = 1/(Δy)^2 .*spdiagm(q, q+2,0 =>ones(T, q),1=> -2 .*ones(T, q), 2 =>ones(T, q))

    xy = Iy'*((Rr ⊗ (∂y ⊗ ∂x)))
    yy = Iy'*((Rr ⊗ (∂y² ⊗ Rp)) + k_b^2 .*Iy)
    zy = Iy'*((∂z ⊗ (∂y ⊗ Rp)))

    #operators for Ez
    Rp = spdiagm(p+1, p+1, 0=>ones(T, p+1))
    Rq = spdiagm(q+1, q+1, 0=>ones(T, q+1))
    Rr = spdiagm(r, r+2, 1=>ones(T, r))

    Iz   = (Rr ⊗ (Rq ⊗ Rp))
    ∂x   = 1/Δx .*spdiagm(p+1, p+2, 0 =>-one(T) .*ones(T, p+1), 1 =>ones(T, p+1))
    ∂y   = 1/Δy .*spdiagm(q+1, q+2, 0 =>-one(T) .*ones(T, q+1), 1 =>ones(T, q+1))
    ∂z   = 1/Δz .*spdiagm(r, r+1, 0 =>-one(T) .*ones(T, r), 1 =>ones(T, r))
    ∂z²  = 1/(Δz)^2 .*spdiagm(r, r+2, 0 =>ones(T, r),1=> -2 .*ones(T, r), 2 =>ones(T, r))

    xz = Iz'*((∂z ⊗ (Rq ⊗ ∂x)))
    yz = Iz'*((∂z ⊗ (∂y ⊗ Rp)))
    zz = Iz'*((∂z² ⊗ (Rq ⊗ Rp)) + k_b^2 .*Iz)

    return hcat(vcat(xx, xy, xz), vcat(yx, yy, yz), vcat(zx, zy, zz));
end

"""
Ig = _create_Greens_functions_restriction_operators(number_of_cells)\n
These operators are used to restrict the output of ∫GχEdV.\n
This integral is computed using an FFT on a domain that is larger the number of cells.\n
Therefore, the output of these should be restricted to the correct Yee cell grid.\n
"""
function _create_Greens_functions_restriction_operators(number_of_cells)
    T = ComplexF32;
    p, q, r = number_of_cells
    dx, dy, dz = 2 .^(ceil.(Int, log2.(number_of_cells.+2).+1))

    Ix = spdiagm(p+2, dx, 0=>ones(T, p+2))
    Iy = spdiagm(q+1, dy, 0=>ones(T, q+1))
    Iz = spdiagm(r+1, dz, 0=>ones(T, r+1))
    x = (Iz ⊗ (Iy ⊗ Ix))

    Ix = spdiagm(p+1, dx, 0=>ones(T, p+1))
    Iy = spdiagm(q+2, dy, 0=>ones(T, q+2))
    y = (Iz ⊗ (Iy ⊗ Ix))

    Iy = spdiagm(q+1, dy, 0=>ones(T, q+1))
    Iz = spdiagm(r+2, dz, 0=>ones(T, r+2))
    z = (Iz ⊗ (Iy ⊗ Ix))

    return blockdiag(x, y, z)
end

"""
Ig = _create_Greens_functions_restriction_operators_gpu(number_of_cells)\n
These operators are used to restrict the output of ∫GχEdV.\n
This integral is computed using an FFT on a domain that is larger the number of 2*cells.\n
Therefore, the output of these should be restricted to the correct Yee cell grid.\n
"""
function _create_Greens_functions_restriction_operators_gpu(number_of_cells)
    T = ComplexF32;
    p, q, r = number_of_cells
    dx, dy, dz = 2 .^(ceil.(Int, log2.(number_of_cells.+2).+1))

    Ix = spdiagm(p+2, dx, 0=>ones(T, p+2))
    Iy = spdiagm(q+1, dy, 0=>ones(T, q+1))
    Iz = spdiagm(r+1, dz, 0=>ones(T, r+1))
    x = (Iz ⊗ (Iy ⊗ Ix))

    Ix = spdiagm(p+1, dx, 0=>ones(T, p+1))
    Iy = spdiagm(q+2, dy, 0=>ones(T, q+2))
    y = (Iz ⊗ (Iy ⊗ Ix))

    Iy = spdiagm(q+1, dy, 0=>ones(T, q+1))
    Iz = spdiagm(r+2, dz, 0=>ones(T, r+2))
    z = (Iz ⊗ (Iy ⊗ Ix))

    Ig = blockdiag(x, y, z) |> CUDA.CUSPARSE.CuSparseMatrixCSC

    return Ig
end

#=
some code not directly required for the method that this package was developed for.
=#
"""
FToE = _magnetic_vector_potential_to_electric_field_operator(number_of_cells, resolution, k_b)\n
This operators is used to compute the total electric field on the FDTD grid (Yee cell) given the magnetic vector potential.\n
k_b is the wavenumber in the background medium.\n
"""
function _magnetic_vector_potential_to_electric_field_operator(
    number_of_cells,
    resolution,
    k_b
)
    T = ComplexF32
    p = number_of_cells[1]
    q = number_of_cells[2]
    r = number_of_cells[3]
    Δx = Float32(resolution[1])
    Δy = Float32(resolution[2])
    Δz = Float32(resolution[3])

    k_b = Float32(k_b)
    # needed -> Hx : (P+1) x (Q+2) : (R+2), wanted -> Ex : P x (Q+1) : (R+1)
    # needed -> Hy : (P+2) x (Q+1) : (R+2), wanted -> Ey : (P+1) x Q : (R+1)
    # needed -> Hz : (P+2) x (Q+2) : (R+1), wanted -> Ez : (P+1) x (Q+1) : R
    # from the new dimension we can map back to the larger E field and get the old solver to work

    #operators for Kx to E

    # an operator to go from the smaller E field domain to the larger one
    Rp = spdiagm(p, p+2, 1 => ones(T, p))
    Rq = spdiagm(q+1, q+1, 0 => ones(T, q+1))
    Rr = spdiagm(r+1, r+1, 0 => ones(T, r+1))

    Ix   = (Rr ⊗ (Rq ⊗ Rp))

    #Kx to E
    ∂y   = 1/Δy .*spdiagm(q+1, q+2, 0 => -one(T) .*ones(T, q+1), 1 => ones(T, q+1))
    ∂z   = 1/Δz .*spdiagm(r+1, r+2, 0 => -one(T) .*ones(T, r+1), 1 => ones(T, r+1))

    xx = spzeros(T, (p+2)*(q+1)*(r+1),(p+1)*(q+2)*(r+2))
    yx = Ix'*((∂z ⊗ (Rq ⊗ Rp))) #f.y -> e.x
    zx = Ix'*((Rr ⊗ (∂y ⊗ Rp))) #f.z -> e.x

    #operators for Ey
    Rp = spdiagm(p+1, p+1, 0 => ones(T, p+1))
    Rq = spdiagm(q, q+2, 1 => ones(T, q))
    Rr = spdiagm(r+1, r+1, 0 => ones(T, r+1))

    Iy   = (Rr ⊗ (Rq ⊗ Rp))

    ∂x   = 1/Δx .*spdiagm(p+1, p+2, 0 => -one(T) .*ones(T, p+1), 1 => ones(T, p+1))
    ∂z   = 1/Δz .*spdiagm(r+1, r+2, 0 => -one(T) .*ones(T, r+1), 1 => ones(T, r+1))

    yy = spzeros(T, (p+1)*(q+2)*(r+1),(p+2)*(q+1)*(r+2))
    xy = Iy'*((∂z ⊗ (Rq ⊗ Rp))) #f.x -> e.y
    zy = Iy'*((Rr ⊗ (Rq ⊗ ∂x))) #f.z -> e.y

    #operators for Ez
    Rp = spdiagm(p+1, p+1, 0 => ones(T, p+1))
    Rq = spdiagm(q+1, q+1, 0=> ones(T, q+1))
    Rr = spdiagm(r, r+2, 1 => ones(T, r))

    Iz   = (Rr ⊗ (Rq ⊗ Rp))

    ∂x   = 1/Δx .*spdiagm(p+1, p+2, 0 => -one(T) .*ones(T, p+1), 1 => ones(T, p+1))
    ∂y   = 1/Δy .*spdiagm(q+1, q+2, 0 => -one(T) .*ones(T, q+1), 1 => ones(T, q+1))

    zz = spzeros(T, (p+1)*(q+1)*(r+2),(p+2)*(q+2)*(r+1))
    xz = Iz'*((Rr ⊗ (∂y ⊗ Rp)))
    yz = Iz'*((Rr ⊗ (Rq ⊗ ∂x)))

    return hcat(vcat(xx, xy, xz), vcat(yx, yy, yz), vcat(zx, zy, zz));
end

"""
AToH = _electric_vector_potential_to_magnetic_field_operator(number_of_cells, resolution, k_b)\n
This operators is used to compute the total magnetic field on the FDTD grid (Yee cell) given the electric vector potential.\n
k_b is the wavenumber in the background medium.\n
"""
function _electric_vector_potential_to_magnetic_field_operator(
    number_of_cells,
    resolution,
    k_b
)
    T = ComplexF32
    p = number_of_cells[1]
    q = number_of_cells[2]
    r = number_of_cells[3]
    Δx = Float32(resolution[1])
    Δy = Float32(resolution[2])
    Δz = Float32(resolution[3])

    k_b = Float32(k_b)
    # needed -> Ex : (P+2) x (Q+1) : (R+1), wanted -> Hx : (P+1) x (Q+2) X (R+2)
    # needed -> Ey : (P+1) x (Q+2) : (R+1), wanted -> Hy : (P+2) x (Q+1) : (R+2)
    # needed -> Ez : (P+1) x (Q+1) : (R+2), wanted -> Hz : (P+2) x (Q+2) : (R+1)
    # from the new dimension we can map back to the larger E field and get the old solver to work

    #operators for Kx to E

    # an operator to go from the smaller E field domain to the larger one
    Rp = spdiagm(p+1, p+1, 0 => ones(T, p+1))
    Rq = spdiagm(q, q+2, 1 => ones(T, q))
    Rr = spdiagm(r, r+2, 1 => ones(T, r))

    Ix   = (Rr ⊗ (Rq ⊗ Rp))

    #Kx to E
    ∂y   = 1/Δy .*spdiagm(q, q+1, 0 => -one(T) .*ones(T, q), 1 => ones(T, q))
    ∂z   = 1/Δz .*spdiagm(r, r+1, 0 => -one(T) .*ones(T, r), 1 => ones(T, r))

    xx = spzeros(T, (p+1)*(q+2)*(r+2),(p+2)*(q+1)*(r+1))
    yx = Ix'*((∂z ⊗ (Rq ⊗ Rp))) #A.y -> h.x
    zx = Ix'*((Rr ⊗ (∂y ⊗ Rp))) #A.z -> h.x

    #operators for Ey
    Rp = spdiagm(p, p+2, 1 => ones(T, p))
    Rq = spdiagm(q+1, q+1, 0 => ones(T, q+1))
    Rr = spdiagm(r, r+2, 1 => ones(T, r))

    Iy   = (Rr ⊗ (Rq ⊗ Rp))

    ∂x   = 1/Δx .*spdiagm(p, p+1, 0 => -one(T) .*ones(T, p), 1 => ones(T, p))
    ∂z   = 1/Δz .*spdiagm(r, r+1, 0 => -one(T) .*ones(T, r), 1 => ones(T, r))

    yy = spzeros(T, (p+2)*(q+1)*(r+2),(p+1)*(q+2)*(r+1))
    xy = Iy'*((∂z ⊗ (Rq ⊗ Rp))) #A.x -> h.y
    zy = Iy'*((Rr ⊗ (Rq ⊗ ∂x))) #A.z -> h.y

    #operators for Ez
    Rp = spdiagm(p, p+2, 1 => ones(T, p))
    Rq = spdiagm(q, q+2, 1 => ones(T, q))
    Rr = spdiagm(r+1, r+1, 0 => ones(T, r+1))

    Iz   = (Rr ⊗ (Rq ⊗ Rp))

    ∂x   = 1/Δx .*spdiagm(p, p+1, 0 => -one(T) .*ones(T, p), 1 => ones(T, p))
    ∂y   = 1/Δy .*spdiagm(q, q+1, 0 => -one(T) .*ones(T, q), 1 => ones(T, q))

    zz = spzeros(T, (p+2)*(q+2)*(r+1),(p+1)*(q+1)*(r+2))
    xz = Iz'*((Rr ⊗ (∂y ⊗ Rp)))
    yz = Iz'*((Rr ⊗ (Rq ⊗ ∂x)))

    return hcat(vcat(xx, xy, xz), vcat(yx, yy, yz), vcat(zx, zy, zz));
end

"""
FToH = _electric_vector_potential_to_magnetic_field_operator(number_of_cells, resolution, k_b)\n
This operators is used to compute the total magnetic field on the FDTD grid (Yee cell) given the magnetic vector potential.\n
k_b is the wavenumber in the background medium.\n
"""
function _magnetic_vector_potential_to_magnetic_field_operator(
    number_of_cells,
    resolution,
    k_b
)
    T = ComplexF32

    p = number_of_cells[1]
    q = number_of_cells[2]
    r = number_of_cells[3]
    Δx = Float32(resolution[1])
    Δy = Float32(resolution[2])
    Δz = Float32(resolution[3])

    k_b = Float32(k_b)
    # needed -> Hx : (P+1) x (Q+2) : (R+2), wanted -> Hx : (P-1) x (Q) : (R)
    # needed -> Hy : (P+2) x (Q+1) : (R+2), wanted -> Hy : (P) x (Q-1) : (R)
    # needed -> Hz : (P+2) x (Q+2) : (R+1), wanted -> Hz : (P) x (Q) : (R-1)
    #operators for Hx
    Rp = spdiagm(p-1, p+1, 1 => ones(T, p-1))
    Rq = spdiagm(q, q+2, 1 => ones(T, q))
    Rr = spdiagm(r, r+2, 1 => ones(T, r))

    Ix   = (Rr ⊗ (Rq ⊗ Rp))
    ∂x   = 1/Δx .*spdiagm(p-1, p+2, 1 => -one(T) .*ones(T, p-1), 2 => ones(T, p-1)) # this could be broken AF ?

    ∂y   = 1/Δy .*spdiagm(q, q+1, 0 => -one(T) .*ones(T, q), 1 => ones(T, q))
    ∂z   = 1/Δz .*spdiagm(r, r+1, 0 => -one(T) .*ones(T, r), 1 => ones(T, r))
    ∂x²  = 1/(Δx)^2 .*spdiagm(p-1, p+1, 0 => ones(T, p-1), 1 => -2 .*ones(T, p-1), 2 => ones(T, p-1))

    xx = Ix'*((Rr ⊗ (Rq ⊗ ∂x²)) + k_b^2 .*Ix) #f.x -> h.x
    yx = Ix'*((Rr ⊗ (∂y ⊗ ∂x))) # f.y -> h.x (P+2) x (Q+1) : (R+2) -> (P-1) x (Q) : (R)
    zx = Ix'*((∂z ⊗ (Rq ⊗ ∂x))) # f.z -> h.x (P+2) x (Q+2) : (R+1) -> (P-1) x (Q) : (R)

    #operators for Ey
    Rp = spdiagm(p, p+2, 1 => ones(T, p))
    Rq = spdiagm(q-1, q+1, 1 => ones(T, q-1))
    Rr = spdiagm(r, r+2, 1 => ones(T, r))

    Iy   = (Rr ⊗ (Rq ⊗ Rp))
    ∂x   = 1/Δx .*spdiagm(p, p+1, 0 => -one(T) .*ones(T, p), 1 => ones(T, p))
    ∂y   = 1/Δy .*spdiagm(q-1, q+2, 1 => -one(T) .*ones(T, q-1), 2 => ones(T, q-1))
    ∂z   = 1/Δz .*spdiagm(r, r+1, 0 => -one(T) .*ones(T, r), 1 => ones(T, r))
    ∂y²  = 1/(Δy)^2 .*spdiagm(q-1, q+1, 0 => ones(T, q-1), 1 => -2 .*ones(T, q-1), 2 => ones(T, q-1))

    xy = Iy'*((Rr ⊗ (∂y ⊗ ∂x)))
    yy = Iy'*((Rr ⊗ (∂y² ⊗ Rp)) + k_b^2 .*Iy)
    zy = Iy'*((∂z ⊗ (∂y ⊗ Rp)))

    #operators for Ez
    Rp = spdiagm(p, p+2, 1 => ones(T, p))
    Rq = spdiagm(q, q+2, 1 => ones(T, q))
    Rr = spdiagm(r-1, r+1, 1 => ones(T, r-1))

    Iz   = (Rr ⊗ (Rq ⊗ Rp))
    ∂x   = 1/Δx .*spdiagm(p, p+1, 0 => -one(T) .*ones(T, p), 1 => ones(T, p))
    ∂y   = 1/Δy .*spdiagm(q, q+1, 0 => -one(T) .*ones(T, q), 1 => ones(T, q))
    ∂z   = 1/Δz .*spdiagm(r-1, r+2, 1 => -one(T) .*ones(T, r-1), 2 => ones(T, r-1))
    ∂z²  = 1/(Δz)^2 .*spdiagm(r-1, r+1, 0 => ones(T, r-1), 1 => -2 .*ones(T, r-1), 2 => ones(T, r-1))

    xz = Iz'*((∂z ⊗ (Rq ⊗ ∂x)))
    yz = Iz'*((∂z ⊗ (∂y ⊗ Rp)))
    zz = Iz'*((∂z² ⊗ (Rq ⊗ Rp)) + k_b^2 .*Iz)

    return hcat(vcat(xx, xy, xz), vcat(yx, yy, yz), vcat(zx, zy, zz));
end

"""
Ig = _create_Greens_functions_restriction_operators_magnetic_field(number_of_cells)\n
These operators are used to restrict the output of ∫GKdV.\n
This integral is computed using an FFT on a domain that is larger the number of 2*cells.\n
Therefore, the output of these should be restricted to the correct Yee cell grid.\n
"""
function _create_Greens_functions_restriction_operators_magnetic_field(number_of_cells)
    T = ComplexF32
    p, q, r = number_of_cells
    dx, dy, dz = 2 .^(ceil.(Int, log2.(number_of_cells.+2).+1))

    Ix = spdiagm(p+1, dx, 0 => ones(T, p+1))
    Iy = spdiagm(q+2, dy, 0 => ones(T, q+2))
    Iz = spdiagm(r+2, dz, 0 => ones(T, r+2))
    x = (Iz ⊗ (Iy ⊗ Ix))

    Ix = spdiagm(p+2, dx, 0 => ones(T, p+2))
    Iy = spdiagm(q+1, dy, 0 => ones(T, q+1))
    y = (Iz ⊗ (Iy ⊗ Ix))

    Iy = spdiagm(q+2, dy, 0 => ones(T, q+2))
    Iz = spdiagm(r+1, dz, 0 => ones(T, r+1))
    z = (Iz ⊗ (Iy ⊗ Ix))

    return blockdiag(x, y, z)
end

"""
Ig = _create_Greens_functions_restriction_operators_magnetic_field_gpu(number_of_cells)\n
These operators are used to restrict the output of ∫GKdV.\n
This integral is computed using an FFT on a domain that is larger the number of 2*cells.\n
Therefore, the output of these should be restricted to the correct Yee cell grid.\n
"""
function _create_Greens_functions_restriction_operators_magnetic_field_gpu(number_of_cells)
    T = ComplexF32
    p, q, r = number_of_cells
    dx, dy, dz = 2 .^(ceil.(Int, log2.(number_of_cells.+2).+1))

    Ix = spdiagm(p+1, dx, 0 => ones(T, p+1))
    Iy = spdiagm(q+2, dy, 0 => ones(T, q+2))
    Iz = spdiagm(r+2, dz, 0 => ones(T, r+2))
    x = (Iz ⊗ (Iy ⊗ Ix))

    Ix = spdiagm(p+2, dx, 0 => ones(T, p+2))
    Iy = spdiagm(q+1, dy, 0 => ones(T, q+1))
    y = (Iz ⊗ (Iy ⊗ Ix))

    Iy = spdiagm(q+2, dy, 0 => ones(T, q+2))
    Iz = spdiagm(r+1, dz, 0 => ones(T, r+1))
    z = (Iz ⊗ (Iy ⊗ Ix))

    Ig = blockdiag(x, y, z) |> CUDA.CUSPARSE.CuSparseMatrixCSC

    return Ig
end


text:
module ImgTransform
    using ImageTransformations
    using ImageBinarization
    using Images

    export ImageTransformations, Images, cutoutShape, imgToBitMatrix

    struct Rectangle
        xl::Int64 # x coordinate of top left corner
        yl::Int64 # y ccordinate of top left corner
        xr::Int64 # x coordinate of bottom right corner
        yr::Int64 # y coordinate of bottom right corner
    end

    function escribedRectangle(binarized::BitMatrix)
        filter = p -> p ≠ 0 # filter used to find relevant points

        h, w = size(binarized)
        rows = [sum(binarized[i, :]) for i in 1:h] # map binarized to the vector of sums of rows, i.e. rows[i] = sum of elements in the i-th row
        columns = [sum(binarized[:, i]) for i in 1:w] # same as before but for columns

        # find relevant points/lines
        top = findfirst(filter, rows) # get index of first non-zero element
        bottom = findlast(filter, rows)
        left = findfirst(filter, columns)
        right = findlast(filter, columns)

        Rectangle(left, top, right, bottom)
    end

    function cutoutShape(binarized::BitMatrix)
        rect = escribedRectangle(binarized)
        binarized[rect.yl:rect.yr, rect.xl:rect.xr] # cut out the rectangle containing the drawing
    end

    function imgToBitMatrix(img::Matrix, negative::Bool = true)
        noalfa = RGB{N0f8}.(img)
        binarized = binarize(noalfa, Otsu()) # binarize using Otsu algorithm; returns Matrix{Gray}
        binarized = BitMatrix(binarized) # convert Matrix{Gray} to BitMatrix
        if negative
            binarized = BitMatrix(xor.(binarized, 1)) # make "negative image"; binarize function transform white color to 1 and black to 0, I want it the other way around
        end

        return binarized
    end

    toVector(m::BitMatrix) = reshape(m, length(m), 1)[1:end]

end

max_stars_repo_name: achjaj/shape-recognition

max_issues_repo_name: achjaj/shape-recognition

max_forks_repo_name: achjaj/shape-recognition