v0.6-Feature Dev

src/Integrals/Engines/Numerical.jl

@@ -1,90 +1,98 @@
 using HCubature
 
-function changeFuncRange(f::F, ::Type{T}) where {F, T}
-    fCore = function (x)
-        x = Tuple(x)
-        f(x ./ (one(T) .- x.^2)) * mapreduce(*, x) do t
-            (1 + t^2) / (1 - t^2)^2
-        end
-    end
-    ReturnTyped(fCore, T)
+struct ConfineInterval{T, F<:ReturnTyped{T}} <: FunctionModifier
+    f::F
+
+    ConfineInterval(f::F) where {T, F<:ReturnTyped{T}} = new{T, F}(f)
 end
 
+ConfineInterval(::Type{T}, f::Function) where {T} = ConfineInterval(ReturnTyped(f, T))
 
-function numericalIntegration(integrand::ReturnTyped{T}, 
-                              interval::NTuple{2, NonEmptyTuple{T, D}}) where {T, D}
-    (first∘hcubature)(integrand, first(interval), last(interval); maxevals=typemax(Int))
+function (f::ConfineInterval{T})(x::AbstractVector{T}) where {T}
+    val = f.f(x ./ (one(T) .- x.^2))
+    mapreduce(*, x) do t
+        tSquare = t * t
+        (1 + tSquare) / (1 - tSquare)^2
+    end * val
 end
 
 
-function composeOneBodyKernel(op::O, ::Type{T}, termL::F1, termR::F2) where 
-                             {O<:Function, T, F1<:Function, F2<:Function}
-    PairCombine(StableBinary(*, T), adjoint∘termL, op(termR))
-end
+struct ModSquaredMap{T, F<:ReturnTyped{T}} <: FunctionModifier
+    f::F
 
-function composeOneBodyKernel(op::O, ::Type{T}, term::F) where {O<:Function, T, F<:Function}
-    composeOneBodyKernel(op, T, term, term)
+    ModSquaredMap(f::F) where {T, F<:ReturnTyped{T}} = new{T, F}(f)
 end
 
-function composeOneBodyKernel(::Identity, ::Type{T}, term::F) where {T, F<:Function}
-    function (input::I) where {I}
-        val = term(input)::T
-        val' * val
-    end
-end
+ModSquaredMap(::Type{T}, f::Function) where {T} = ModSquaredMap(ReturnTyped(f, T))
 
-function numericalOneBodyInt(op::F, (orb,)::Tuple{EvalFieldFunction{T, D}}, 
-                             (pVal,)::Tuple{FilteredVecOfArr{T}}) where 
-                            {F<:DirectOperator, T, D}
-    term = Base.Fix2(orb, pVal)
-    integrand = changeFuncRange(composeOneBodyKernel(op, T, term), T)
-    bound = ntuple(_->one(T), Val(D))
-    numericalIntegration(integrand, (.-(bound), bound))
+function (f::ModSquaredMap{T})(input::AbstractVector{T}) where {T}
+    val = f.f(input)
+    val' * val
 end
 
-function numericalOneBodyInt(op::F, orbs::NTuple{2, EvalFieldFunction{T, D}}, 
-                             pValPair::NTuple{2, FilteredVecOfArr{T}}) where 
-                            {F<:DirectOperator, T, D}
-    termL, termR = Base.Fix2.(orbs, pValPair)
-    integrand = changeFuncRange(composeOneBodyKernel(op, T, termL, termR), T)
-    bound = ntuple(_->one(T), Val(D))
-    numericalIntegration(integrand, (.-(bound), bound))
+
+function numericalIntegrateCore(integrand::ConfineInterval{T}, 
+                                interval::NTuple{2, NonEmptyTuple{T, D}}) where {T, D}
+    fullRes = hcubature(integrand, first(interval), last(interval), maxevals=typemax(Int))
+    first(fullRes)
 end
 
 
-struct NumOverlapOrbSelf{T, D, B<:PrimitiveOrbCore{T, D}} <: OrbitalIntegrator{T, D}
-    orb::Tuple{B}
+function composeOneBodyKernel(op::O, ::Type{T}, (termL, termR)::Tuple{F1, F2}) where 
+                             {O<:Function, T, F1<:Function, F2<:Function}
+    PairCombine(StableMul(T), adjoint∘termL, op(termR))
 end
 
-function (f::NumOverlapOrbSelf{T})(pars::FilteredVecOfArr{T}) where {T}
-    numericalOneBodyInt(Identity(), f.orb, (pars,))
+function composeOneBodyKernel(op::O, ::Type{T}, (term,)::Tuple{F}) where 
+                             {O<:Function, T, F<:Function}
+    composeOneBodyKernel(op, T, (term, term))
 end
 
-struct NumOverlapOrbPair{T, D, B1<:PrimitiveOrbCore{T, D}, 
-                         B2<:PrimitiveOrbCore{T, D}} <: OrbitalIntegrator{T, D}
-    orb::Tuple{B1, B2}
+function composeOneBodyKernel(::Identity, ::Type{T}, (term,)::Tuple{F}) where 
+                             {T, F<:Function}
+    ModSquaredMap(T, term)
 end
 
-function (f::NumOverlapOrbPair{T})(pars::Vararg{FilteredVecOfArr{T}, 2}) where {T}
-    numericalOneBodyInt(Identity(), f.orb, pars)
+
+function composeIntegralKernel(::OneBodyIntegral, op::O, ::Type{T}, 
+                               terms::N12Tuple{Function}) where {O<:Function, T}
+    composeOneBodyKernel(op, T, terms)
 end
 
-genOneBodyCoreIntegrator(::Identity, orbs::Tuple{PrimitiveOrbCore{T, D}}) where {T, D} = 
-NumOverlapOrbSelf(orbs)
+function numericalIntegrate(::OneBodyIntegral{D}, op::F, 
+                            orbs::NonEmptyTuple{EvalFieldFunction{T, D}, N}, 
+                            pVals::NonEmptyTuple{FilteredVecOfArr{T}, N}) where 
+                            {F<:DirectOperator, T, D, N}
+    terms = Base.Fix2.(orbs, pVals)
+    integralKernel = composeIntegralKernel(OneBodyIntegral{D}(), op, T, terms)
+    integrand = ConfineInterval(T, integralKernel)
+    bound = ntuple(_->one(T), Val(D))
+    numericalIntegrateCore(integrand, (.-(bound), bound))
+end
 
-genOneBodyCoreIntegrator(::Identity, orbs::NTuple{2, PrimitiveOrbCore{T, D}}) where {T, D} = 
-NumOverlapOrbPair(orbs)
 
 function buildNormalizerCore(o::PrimitiveOrbCore{T, D}) where {T, D}
-    genOneBodyCoreIntegrator(Identity(), (o,))
+    buildCoreIntegrator(OneBodyIntegral{D}(), Identity(), (o,))
+end
+
+
+struct OneBodyNumIntegrate{T, D, F<:DirectOperator, 
+                           P<:N12Tuple{PrimitiveOrbCore{T, D}}} <: OrbitalIntegrator{T, D}
+    op::F
+    basis::P
 end
 
+const OneBodySelfNumInt{T, D, F<:DirectOperator, P<:PrimitiveOrbCore{T, D}} = 
+      OneBodyNumIntegrate{T, D, F, Tuple{P}}
+const OnyBodyPairNumInt{T, D, F<:DirectOperator, P1<:PrimitiveOrbCore{T, D}, 
+                        P2<:PrimitiveOrbCore{T, D}} = 
+      OneBodyNumIntegrate{T, D, F, Tuple{P1, P2}}
+
+function (f::OneBodySelfNumInt{T, D})(pVal::FilteredVecOfArr{T}) where {T, D}
+    numericalIntegrate(OneBodyIntegral{D}(), f.op, f.basis, (pVal,))
+end
 
-#!! Implement a functional struct (<:OrbitalIntegrator{T, D}) for the closure
-function genOneBodyCoreIntegrator(op::F, 
-                                  orbs::NonEmptyTuple{PrimitiveOrbCore{T, D}, N}) where 
-                                 {F<:DirectOperator, N, T, D}
-    function (pVal::AbtVecOfAbtArr{T}, pVals::Vararg{AbtVecOfAbtArr{T}, N}) where {T}
-        numericalOneBodyInt(op, orbs, (pVal, pVals...))
-    end
+function (f::OnyBodyPairNumInt{T, D})(pVal1::FilteredVecOfArr{T}, 
+                                      pVal2::FilteredVecOfArr{T}) where {T, D}
+    numericalIntegrate(OneBodyIntegral{D}(), f.op, f.basis, (pVal1, pVal2))
 end