From 3209078dc5999a1dc5d19750825babea99979427 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Tue, 29 Aug 2023 18:45:45 -0700 Subject: [PATCH 01/22] Developing new parameter containers. --- Project.toml | 2 + src/AbstractTypes.jl | 46 +++- src/Basis.jl | 41 ++-- src/CI.jl | 129 ++++++++++++ src/Exception.jl | 16 ++ src/Nuclear.jl | 26 +++ src/Overload.jl | 27 ++- src/Parameters.jl | 325 ++++++++++++++++++++--------- src/Quiqbox.jl | 1 + src/Tools.jl | 24 ++- test/AngularMomentum-test.jl | 62 ++++++ test/runtests.jl | 2 +- test/unit-tests/Basis-test.jl | 2 +- test/unit-tests/Parameters-test.jl | 2 +- test/unit-tests/Tools-test.jl | 2 +- 15 files changed, 566 insertions(+), 141 deletions(-) create mode 100644 src/CI.jl create mode 100644 src/Exception.jl create mode 100644 src/Nuclear.jl create mode 100644 test/AngularMomentum-test.jl diff --git a/Project.toml b/Project.toml index 24fb0cef..0e56cdd2 100644 --- a/Project.toml +++ b/Project.toml @@ -11,9 +11,11 @@ LBFGSB = "5be7bae1-8223-5378-bac3-9e7378a2f6e6" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" SPGBox = "bf97046b-3e66-4aa0-9aed-26efb7fac769" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" [compat] diff --git a/src/AbstractTypes.jl b/src/AbstractTypes.jl index d2a6a2c3..df948e83 100644 --- a/src/AbstractTypes.jl +++ b/src/AbstractTypes.jl @@ -39,8 +39,10 @@ abstract type ManyFermionDataBox{T, D} <: ImmutableDataBox{T} end abstract type HartreeFockintermediateData{T} <: MutableDataBox{T} end -abstract type DifferentiableParameter{DataT, ContainerT} <: SemiMutableParameter{DataT, ContainerT} end -abstract type ConfigBox{T, ContainerT, MethodT} <: MutableParameter{ContainerT, T} end +abstract type VariableBox{T, ContainerT} <: SemiMutableParameter{T, ContainerT} end +abstract type ConfigBox{T, ContainerT, MethodT} <: MutableParameter{T, ContainerT} end + +abstract type ParamBox{T, V, F} <: VariableBox{T, ParamBox} end abstract type HartreeFockFinalValue{T, HFT} <: AbstractHartreeFockFinalValue{T} end @@ -101,4 +103,42 @@ const NTupleOfSBN{BN, T, D} = NTuple{BN, SpatialBasis{T, D}} const StrOrSym = Union{String, Symbol} -const MatrixCol{T} = Vector{SubArray{T, 1, Matrix{T}, Tuple{Slice{OneTo{Int}}, Int}, true}} \ No newline at end of file +const MatrixCol{T} = Vector{SubArray{T, 1, Matrix{T}, Tuple{Slice{OneTo{Int}}, Int}, true}} + +const Array0D{T} = Array{T, 0} + +const Union0D{T} = Union{T, Array0D{T}} + +const MTuple{M, T} = NTuple{M, Array0D{T}} # M: mutable + +const UTuple{U, T} = NTuple{U, Union0D{T}} # U: union + +const NMOTuple{NMO, T} = Tuple{T, Vararg{T, NMO}} + +const MMOTuple{MMO, T} = Tuple{Array0D{T}, Vararg{Array0D{T}, MMO}} + +const UMOTuple{UMO, T} = Tuple{Union0D{T}, Vararg{Union0D{T}, UMO}} + +const NMTTuple{NMO, T} = Tuple{T, T, Vararg{T, NMO}} + +const MMTTuple{MMO, T} = Tuple{Array0D{T}, Array0D{T}, Vararg{Array0D{T}, MMO}} + +const UMTTuple{UMO, T} = Tuple{Union0D{T}, Union0D{T}, Vararg{Union0D{T}, UMO}} + +const TorTuple{T} = Union{T, Tuple{T}} + +const TorTupleLong{T} = Union{T, Tuple{T, Vararg{T}}} + +const TorAbstractVec{T} = Union{T, AbstractVector{<:T}} + +const IntOrNone = Union{Int, Nothing} +const IntOrMiss = Union{Int, Missing} + +const MutableIndex = Array0D{IntOrNone} + + +const MonoParam{T} = OneParam{T} + +const DualParam{T} = Tuple{OneParam{T}, OneParam{T}} + +const MoreParam{T} = Tuple{OneParam{T}, OneParam{T}, OneParam{T}, Vector{OneParam{T}}} \ No newline at end of file diff --git a/src/Basis.jl b/src/Basis.jl index 8b78c6e3..61fa1e44 100644 --- a/src/Basis.jl +++ b/src/Basis.jl @@ -68,24 +68,25 @@ getTypeParams(::GaussFunc{T, Fxpn, Fcon}) where {T, Fxpn, Fcon} = (T, Fxpn, Fcon {T<:AbstractFloat} -> ParamBox{T, :$(xpnSym)} -Construct an exponent coefficient given a value or variable (with its symbol). -`mapFunction`, `canDiff`, and `inSym` work the same way as in a general constructor of a -[`ParamBox`](@ref). +Construct an exponent coefficient given a value or variable (and its symbol without the +index). `mapFunction`, `canDiff`, and `inSym` work the same way as in a general constructor +of a [`ParamBox`](@ref). """ -genExponent(eData::Pair{Array{T, 0}, Symbol}, mapFunction::Function=itself; +genExponent(eData::TorTuple{OneParam{T}}, mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true)) where {T<:AbstractFloat} = -ParamBox(Val(xpnSym), mapFunction, eData, genIndex(nothing), fill(canDiff)) +ParamBox(Val(xpnSym), mapFunction, eData, genIndex(), fill(canDiff)) -genExponent(e::Array{T, 0}, mapFunction::Function=itself; +genExponent(eData::TorTuple{Union0D{T}}, mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true), - inSym::Symbol=xpnIVsym) where {T<:AbstractFloat} = -genExponent(e=>inSym, mapFunction; canDiff) + inSym::TorTuple{Symbol}=(xpnIVsym,)) where {T<:AbstractFloat} = +genExponent(OneParam.(tupleObj(eData), inSym), mapFunction; canDiff) -genExponent(e::AbstractFloat, mapFunction::Function=itself; +genExponent(eData::NMOTuple{NPMO, Union{Union0D{T}, OneParam{T}}}, + mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true), - inSym::Symbol=xpnIVsym) = -genExponent(fill(e)=>inSym, mapFunction; canDiff) + inSym::Symbol=xpnIVsym) where {T<:AbstractFloat} = +genExponent(OneParam.(tupleObj(eData), inSym), mapFunction; canDiff) """ @@ -114,15 +115,15 @@ genExponent(pb::ParamBox) = ParamBox(Val(xpnSym), pb) {T<:AbstractFloat} -> ParamBox{T, :$(conSym)} -Construct a contraction coefficient given a value or variable (with its symbol). -`mapFunction`, `canDiff`, and `inSym` work the same way as in a general constructor of a -[`ParamBox`](@ref). +Construct a contraction coefficient given a value or variable (and its symbol without the +index).`mapFunction`, `canDiff`, and `inSym` work the same way as in a general constructor +of a [`ParamBox`](@ref). """ genContraction(dData::Pair{Array{T, 0}, Symbol}, mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true)) where {T<:AbstractFloat} = -ParamBox(Val(conSym), mapFunction, dData, genIndex(nothing), fill(canDiff)) +ParamBox(Val(conSym), mapFunction, dData, genIndex(), fill(canDiff)) genContraction(d::Array{T, 0}, mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true), @@ -268,7 +269,7 @@ genSpatialPoint(compData::Pair{Array{T, 0}, Symbol}, compIndex::Int, mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true)) where {T<:AbstractFloat} = -ParamBox(Val(SpatialParamSyms[compIndex]), mapFunction, compData, genIndex(nothing), +ParamBox(Val(SpatialParamSyms[compIndex]), mapFunction, compData, genIndex(), fill(canDiff)) genSpatialPoint(comp::Array{T, 0}, compIndex::Int, @@ -1868,9 +1869,9 @@ end Return the parameter(s) stored in the input container. If `symbol` is set to `missing`, then return all the parameter(s). If it's set to a `Symbol` tied to a parameter, for example `:α₁`, the function will match any `pb::`[`ParamBox`](@ref) such that -[`indVarOf`](@ref)`(pb)[begin] == :α₁`. If it's set to a `Symbol` without any subscript, +[`indSymOf`](@ref)`(pb) == :α₁`. If it's set to a `Symbol` without any subscript, for example `:α`, the function will match it with all the `pb`s such that -`string(indVarOf(pb)[begin])` contains `'α'`. If the first argument is a collection, its +`string(indSymOf(pb))` contains `'α'`. If the first argument is a collection, its entries must be `ParamBox` containers. """ getParams(pb::ParamBox, symbol::Union{Symbol, Missing}=missing) = @@ -1908,7 +1909,7 @@ getParams(@nospecialize(cs::Tuple{QuiqboxContainer, Vararg{QuiqboxContainer}}), getParams(lazyCollect(cs), symbol) paramFilter(pb::ParamBox, sym::Union{Symbol, Missing}) = -sym isa Missing || inSymbol(sym, indVarOf(pb)[begin]) +sym isa Missing || inSymbol(sym, indSymOf(pb)) copyParContainer(pb::ParamBox, @@ -1989,7 +1990,7 @@ independent variable (in other words, considered identical in the differentiatio procedure) will be marked with same index. `filterMapping` determines whether filtering out (i.e. not return) the extra `ParamBox`es that hold the same independent variables but represent different output variables. The independent variable held by a `ParamBox` can be -inspected using [`indVarOf`](@ref). +inspected using [`indSymOf`](@ref). """ markParams!(b::Union{AbstractVector{T}, T}, filterMapping::Bool=false) where {T<:ParameterizedContainer} = diff --git a/src/CI.jl b/src/CI.jl new file mode 100644 index 00000000..5bfb4c90 --- /dev/null +++ b/src/CI.jl @@ -0,0 +1,129 @@ +struct RSD{𝑚ₛ, ON} <:RCI{𝑚ₛ, ON} + N::Int + spinUpOccu::NTuple{ON, Bool} + spinDnOccu::NTuple{ON, Bool} + ID::UInt + + function RSD(spinUpOccu::NTuple{ON, Bool}, spinDnOccu::NTuple{ON, Bool}, + basis::NTuple{ON, <:GTBasisFuncs{T, D, 1}}) + Nup = sum(spinUpOccu) + Ndn = sum(spinDnOccu) + new{Nup-Ndn, ON}(Nup+Ndn, spinUpOccu, spinDnOccu, objectid(basis)) + end +end + +struct CSF{𝑠, 𝑚ₛ, ON} <: RCI{𝑚ₛ, ON} + N::Int + spatialOccu::NTuple{ON, Int} + ID::UInt + + function RSD(spatialOccu::NTuple{ON, Bool}, basis::NTuple{ON, <:GTBasisFuncs{T, D, 1}}) + + for + new{Nup-Ndn, ON}(Nup+Ndn, spinUpOccu, spinDnOccu, objectid(basis)) + end +end + + + + +function showbitStr(num, nDigits=4) + str = num|>bitstring + str[end-nDigits+1:end] +end + +function computeSpinPerm(nUp, nDn) # potential improvement: only iterate over first half. + u = (1 << nUp) - 1 + siz = binomial(nUp+nDn, nUp) + v = Array{Int}(undef, siz) + v[begin] = u + for i in 2:siz + t = u | (u - 1) + t1 = t + 1 + t2 = ((~t & t1) - 1) >> (trailing_zeros(u) + 1) + u = t1 | t2 + # @show u + # @show showbitStr.([t, t1, t2], nUp+nDn) + v[i] = u + end + # @show showbitStr.(v, nUp+nDn) + v +end + +## Input: [1, 0, 1, 2] +function genSpinStrings(spatialOrbConfig::Vector{Int}) + ids1 = findall(isone, spatialOrbConfig) + ids2 = findall(isequal(2), spatialOrbConfig) + nSingleOccu = length(ids1) + # nDoubleOccu = sum(ids2) + nUp = nSingleOccu ÷ 2 + nDn = nSingleOccu - nUp + permStrs = showbitStr.(computeSpinPerm(nUp, nDn), nSingleOccu) + spinUpConfigs = collect(zero(spatialOrbConfig) for _ in eachindex(permStrs)) + spinDnConfigs = deepcopy(spinUpConfigs) + for (perm, SUconfig, SDconfig) in zip(permStrs, spinUpConfigs, spinDnConfigs) + SUconfig[ids2] .= 1 + SDconfig[ids2] .= 1 + for (spin, idx1) in zip(perm, ids1) + config = spin=='1' ? SUconfig : SDconfig + config[idx1] = 1 + end + end + spinUpConfigs, spinDnConfigs +end + + + +# RHF state +function formatSpatialOrbConfig((RHFSconfig,)::Tuple{NTuple{ON, String}}) + res = (Int[], Int[], Int[]) # [DoubleOccuIds, SingleOccuIds, UnOccuIds] + for (i, e) in enumerate(RHFSconfig) + if e == spinOccupations[end] + push!(res[1], i) + elseif e == spinOccupations[begin] + push!(res[3], i) + else + push!(res[2], i) + end + end + res +end + +function genSpatialOrbConfigCore(n::Int, refConfig::NTuple{3, Vector{Int}}) + N = 2*length(refConfig[begin]) + length(refConfig[begin+1]) + for i in min(n, N) + + end +end + +function cc() + +end + +function promoteElec!(spoConifg, aᵢ, cⱼ) + spoConifg[aᵢ] -= 1 + spoConifg[cⱼ] += 1 + spoConifg +end + +[2,2,2,0,0,0] + + + + +function restrictedSDconfig() + getEhf() + sum(diag(Hc)[1:n]) + 0.5*sum([(eeI[i,i,j,j] - eeI[i,j,j,i]) for j in 1:n, i in 1:n]) +end + +function genFCI() + +end + +function genXCI() + +end + +function genSCI() + +end \ No newline at end of file diff --git a/src/Exception.jl b/src/Exception.jl new file mode 100644 index 00000000..a5823f84 --- /dev/null +++ b/src/Exception.jl @@ -0,0 +1,16 @@ +function checkFuncReturn(f::F, fName::Symbol, ::Type{TArg}, ::Type{TReturn}) where + {F<:Function, TArg, TReturn} + Base.return_types(f, (TArg,))[] == TReturn || + throw(ArgumentError("`$fName`: `$f` should return a `$TReturn` given an input "* + "argument of type `$TArg`.")) +end + + +function checkCollectionMinLen(data, dataSym::Symbol, minLen::Int) + minLen = abs(minLen) + str = minlen > 1 ? "s." : "." + dataLen = length(data) + dataLen < minLen && + throw(ArgumentError("`$dataSym`: $data should contain at least $minLen element"* str)) + dataLen == minLen +end \ No newline at end of file diff --git a/src/Nuclear.jl b/src/Nuclear.jl new file mode 100644 index 00000000..0441ec55 --- /dev/null +++ b/src/Nuclear.jl @@ -0,0 +1,26 @@ +struct Nuc{T, D, NN, CPT} + sym::NTuple{NN, Symbol} + pos::Tuple{Vararg{SpatialPoint{T, D, <:CPT}, NN}} + + function Nuc{CPT}(@nospecialize(sps::Tuple{SpatialPoint{T, D, <:CPT}, + Vararg{SpatialPoint{T, D, <:CPT}, + NN}})) where {CPT, T, D, NN} + sym = map(sps) do sp + marker = sp.marker + haskey(AtomicNumberList, marker) || (throw∘KeyError)(marker) + marker + end + new{T, D, NN+1, CPT}(sym, sps) + end +end + +Nuc(sps::Tuple{SpatialPoint{T, D, CPT}, Vararg{SpatialPoint{T, D, CPT}}}) where + {T, D, CPT} = +Nuc{CPT}(sps) + +Nuc(@nospecialize(sps::Tuple{SpatialPoint{T, D, <:NTuple{D, ParamBox{T}}}, + Vararg{SpatialPoint{T, D, <:NTuple{D, ParamBox{T}}}}})) where + {T, D} = +Nuc{NTuple{D, ParamBox{T}}}(sps) + +Nuc(sps::AbstractVector{<:SpatialPoint{T, D}}) where {T, D} = (Nuc∘Tuple)(sps) \ No newline at end of file diff --git a/src/Overload.jl b/src/Overload.jl index 4d7c8ac1..dddb097c 100644 --- a/src/Overload.jl +++ b/src/Overload.jl @@ -1,9 +1,13 @@ # Julia internal methods overload. import Base: == +==(op1::OneParam, op2::OneParam) = op1.data == op2.data && + op1.label == op2.label && + op1.index[] === op2.index[] + ==(pb1::ParamBox, pb2::ParamBox) = (pb1.data == pb2.data && isDiffParam(pb1) == isDiffParam(pb2) && - pb1.index[] == pb2.index[] && + pb1.index[] === pb2.index[] && pb1.map == pb2.map) ==(cp1::SpatialPoint, cp2::SpatialPoint) = (cp1.param == cp2.param) @@ -85,7 +89,7 @@ end function printIndVarCore(io::IO, pb::ParamBox) - printstyled(io, "$(indVarOf(pb)[begin])", color=:cyan) + printstyled(io, "$(indSymOf(pb))", color=:cyan) end function printIndVar(io::IO, @nospecialize(pb::ParamBox)) @@ -300,12 +304,14 @@ import Base: * ## Iteration Interface import Base: iterate, size, length, eltype -iterate(pb::ParamBox) = (pb.data[][begin][], nothing) -iterate(@nospecialize(_::ParamBox), _) = nothing -size(::ParamBox) = () -length(::ParamBox) = 1 + +iterate(pb::ParamBox) = iterate(pb.data[]) +iterate(pb::ParamBox, i) = iterate(pb.data[], i) +size(::ParamBox{<:Any, <:Any, <:Any, 1}) = () +size(::ParamBox{<:Any, <:Any, <:Any, NP}) where {NP} = (NP,) +size(pb::ParamBox, d::Integer) = size(pb.data[], d) +length(::ParamBox{<:Any, <:Any, <:Any, NP}) where {NP} = NP eltype(::ParamBox{T}) where {T} = T -size(pb::ParamBox, d::Integer) = size(pb.data[][begin][], d) iterate(sp::SpatialPoint) = iterate(sp.param) iterate(sp::SpatialPoint, state) = iterate(sp.param, state) @@ -350,10 +356,13 @@ end ## Indexing Interface import Base: getindex, setindex!, firstindex, lastindex, eachindex, axes -getindex(pb::ParamBox) = pb.data[][begin][] +getindex(op::OneParam) = op.data[] +getindex(pb::ParamBox) = inValOf(pb) getindex(pb::ParamBox, ::Val{:first}) = getindex(pb) getindex(pb::ParamBox, ::Val{:last}) = getindex(pb) -setindex!(pb::ParamBox, d) = begin pb.data[][begin][] = d end +setindex!(op::OneParam, d) = begin op.data[] = d end +setindex!(pb::ParamBox{<:Any, <:Any, <:Any, 1}, d) = begin pb.data[][begin] = d end +setindex!(pb::ParamBox, d, index) = begin pb.data[][index] = d end firstindex(@nospecialize(_::ParamBox)) = Val(:first) lastindex(@nospecialize(_::ParamBox)) = Val(:last) axes(@nospecialize(_::ParamBox)) = () diff --git a/src/Parameters.jl b/src/Parameters.jl index 32b03638..b813558e 100644 --- a/src/Parameters.jl +++ b/src/Parameters.jl @@ -1,6 +1,33 @@ export ParamBox, inValOf, outValOf, inSymOf, outSymOf, isInSymEqual, isOutSymEqual, - indVarOf, dataOf, mapOf, outValCopy, fullVarCopy, enableDiff!, disableDiff!, - isDiffParam, toggleDiff!, changeMapping + indParOf, indSymOf, dataOf, mapOf, outValCopy, fullVarCopy, enableDiff!, + disableDiff!, isDiffParam, toggleDiff!, changeMapping + +using StructArrays + +struct OneParam{T} <: VariableBox{T, OneParam} + label::Symbol + data::Array0D{T} + index::MutableIndex +end + +OneParam(label::Symbol, data::Union0D{T}, index::IntOrNone=nothing) where {T} = +OneParam{T}(label, fillObj(data), genIndex(index)) + +OneParam(labelToData::Pair{Symbol, <:Union0D{T}}, index::IntOrNone=nothing) where {T} = +OneParam{T}(labelToData[begin], fillObj(labelToData[end]), genIndex(index)) + +OneParam(op::OneParam) = itself(op) + +const OPessential{T} = Union{OneParam{T}, Pair{Symbol, <:Union0D{T}}} + +const OneParamStructVec{T} = + StructVector{OneParam{T}, NamedTuple{(:label, :data, :index), + Tuple{Vector{Symbol}, Vector{Array0D{T}}, Vector{MutableIndex}}}, Int} + + +valOf(par::OneParam) = par.data[] +symOf(par::OneParam) = Symbol(par.label, numToSubs(par.index[])) +# indexOf(par::OneParam) = par.index[] """ @@ -48,7 +75,8 @@ stored. If the latter is the type of `data`, then it will directly used to const `outSym::Symbol`: The symbol of the output variable represented by the constructed `ParamBox`. It's equal to the type parameter `V` of the constructed `ParamBox`. -`inSym::Symbol`: The symbol of the input variable held by the constructed `ParamBox`. +`inSym::Symbol`: The symbol of the input variable (without the index) held by the +constructed `ParamBox`. `mapFunction::Function`: The mapping (`mapFunction(::T)->T`) of the input variable, which will be assigned to the field `.map`. When `mapFunction` is not provided, `.map` is set to @@ -88,49 +116,143 @@ parameter functional (e.g., the Hartree–Fock energy). However, the derivative to the stored input variable can also be computed to when the `ParamBox` is marked as differentiable. """ -struct ParamBox{T, V, F<:Function} <: DifferentiableParameter{T, ParamBox} - data::Array{Pair{Array{T, 0}, Symbol}, 0} +struct MonoParamBox{T, V, F<:Function} <: ParamBox{T, V, F} + data::Array0D{OneParam{T}} + map::Union{F, DI{F}} + index::MutableIndex + canDiff::Array0D{Bool} + + function MonoParamBox{T, V}(mapFunction::F, data::OneParam{T}, + index::Array0D{<:IntOrNone}, canDiff::Array0D{Bool}) where + {T, V, F<:Function} + checkFuncReturn(mapFunction, :mapFunction, T, T) + new{T, V, dressOf(F)}(fill(data), mapFunction, index, canDiff) + end + + MonoParamBox{T, V}(::IF, data::OneParam{T}, index, canDiff) where {T, V} = + new{T, V, iT}(fill(data), itself, index, canDiff) +end + +const DefaultMPBoxMap = itself + +struct PolyParamBox{T, V, F<:Function} <: ParamBox{T, V, F} + data1::Array0D{OneParam{T}} + data2::Array0D{OneParam{T}} + data3::OneParamStructVec{T} map::Union{F, DI{F}} - canDiff::Array{Bool, 0} - index::Array{Union{Int, Nothing}, 0} - - function ParamBox{T, V}(f::F, data::Pair{Array{T, 0}, Symbol}, index, canDiff) where - {T, V, F<:Function} - Base.return_types(f, (T,))[1] == T || - throw(AssertionError("The mapping function `f`: `$(f)` should return the same "* - "data type as its input argument.")) - new{T, V, dressOf(F)}(fill(data), f, canDiff, index) + index::MutableIndex + canDiff::Array0D{Bool} + + function PolyParamBox{T, V}(mapFunction::F, data1::OneParam{T}, data2::OneParam{T}, + data3::AbstractVector{OneParam{T}}, + index::Array0D{<:IntOrNone}, canDiff::Array0D{Bool}) where + {T, V, F<:Function} + checkFuncReturn(mapFunction, :mapFunction, Vector{T}, T) + new{T, V, dressOf(F)}(fill(data1), fill(data2) StructArray(data3), + f, index, canDiff) end +end - ParamBox{T, V}(data::Pair{Array{T, 0}, Symbol}, index, canDiff) where {T, V} = - new{T, V, iT}(fill(data), itself, canDiff, index) +const DefaultPPBoxMap = sum + +genParamBox(::Val{V}, mapFunction::F, data::OneParam{T}, + index::Union0D{<:IntOrNone}=genIndex(), + canDiff::Union0D{Bool}=getFLevel(F)>0) where {T, V, F<:Function} = +MonoParamBox{T, V}(mapFunction, data, fillObj(index), fillObj(canDiff)) + +genParamBox(::Val{V}, mapFunction::F, data1::OneParam{T}, data2::OneParam{T}, + data3::AbstractVector{OneParam{T}}, + index::Union0D{<:IntOrNone}=genIndex(), + canDiff::Union0D{Bool}=false) where {T, V, F<:Function} = +PolyParamBox{T, V}(mapFunction, data1, data2, data3, fillObj(index), fillObj(canDiff)) + +function genParamBox(::Val{V}, mapFunction::F, data::AbstractVector{OneParam{T}}, + index::Union0D{<:IntOrNone}=genIndex(), + canDiff::Union0D{Bool}=getFLevel(F)>0) where {T, V, F<:Function} + isLenOne = checkCollectionMinLen(data, :data, 1) + if isLenOne + genParamBox(Val(V), mapFunction, data[], index, canDiff) + else + genParamBox(Val(V), mapFunction, data[begin], data[begin+1], data[begin+2:end], + index, canDiff) + end end -ParamBox(::Val{V}, mapFunction::F, data::Pair{Array{T, 0}, Symbol}, - index=genIndex(nothing), canDiff=fill(true)) where {V, F<:Function, T} = -ParamBox{T, V}(mapFunction, data, index, canDiff) +genParamBox(::Val{V}, pb::MonoParamBox{T}, index::Union0D{<:IntOrNone}=genIndex(), + canDiff::Array0D{Bool}=copy(pb.canDiff)) where {V, T} = +genParamBox(Val(V), pb.map, pb.data[], index, canDiff) + +genParamBox(::Val{V}, pb::PolyParamBox{T}, index::Union0D{<:IntOrNone}=genIndex(), + canDiff::Array0D{Bool}=copy(pb.canDiff)) where {V, T} = +genParamBox(Val(V), pb.map, pb.data1[], pb.data2[], pb.data3, index, canDiff) + +genParamBox(data::OPessential{T}, outSym::Symbol=:undef, mapFunction::F=DefaultMPBoxMap; + index::IntOrNone=nothing, + canDiff::Bool=getFLevel(F)>0) where {T, F<:Function} = +genParamBox(Val(outSym), mapFunction, OneParam(data), index, fill(canDiff)) + +function genParamBox(data::AbstractVector{<:Union{Union0D{T}, OPessential{T}}}, + outSym::Symbol=:undef, + mapFunction::F=ifelse(checkCollectionMinLen(data, :data, 1), + DefaultMPBoxMap, DefaultPPBoxMap); + defaultDataLabel::Symbol=Symbol(IVsymSuffix, outSym), + index::IntOrNone=nothing, + canDiff::Bool=getFLevel(F)>0) where {T, F<:Function} + data = map(data) do item + item isa Union0D ? OneParam(defaultDataLabel, item) : OneParam(item) + end + genParamBox(Val(outSym), mapFunction, data, index, canDiff) +end -ParamBox(::Val{V}, ::IF, data::Pair{Array{T, 0}, Symbol}, index=genIndex(nothing), - canDiff=fill(false)) where {V, T} = -ParamBox{T, V}(data, index, canDiff) +# function ParamBox(data::Union0D{T}, outSym::Symbol=:undef, mapFunction::F=itself; +# defaultDataLabel::Symbol=Symbol(IVsymSuffix, outSym), +# index::IntOrNone=nothing, canDiff::Bool=getFLevel(F)>0) where +# {T, F<:Function} +# data = map(data|>tupleObj) do item +# item isa Union0D ? OneParam(defaultDataLabel, item) : item +# end +# ParamBox(Val(outSym), mapFunction, OneParam.(tupleObj(data)), +# genIndex(index), fill(canDiff)) +# end + +# function ParamBox(data::TorTupleLong{Union{Union0D{T}, OPessential{T}}}, +# outSym::Symbol=:undef, mapFunction::Function=itself; +# defaultDataLabel::Symbol=Symbol(IVsymSuffix, outSym), +# index::Union{Int, Nothing}=nothing, canDiff::Bool=false) +# data = map(data|>tupleObj) do item +# item isa Union0D ? OneParam(defaultDataLabel, item) : item +# end +# ParamBox(Val(outSym), mapFunction, OneParam.(tupleObj(data)), +# genIndex(index), fill(canDiff)) +# end + +# ParamBox(inVar::Union0D{T}, outSym::Symbol=:undef, +# inSym::Symbol=(Symbol(IVsymSuffix, outSym),); +# index::Union{Int, Nothing}=nothing, canDiff::Bool=false) where {T} = +# ParamBox(Val(outSym), itself, OneParam(inVar, inSym), genIndex(index), fill(canDiff)) + +# ParamBox(inVar::Union0D{T}, outSym::Symbol, mapFunction::Function, +# inSym::Symbol=Symbol(IVsymSuffix, outSym); +# index::Union{Int, Nothing}=nothing, canDiff::Bool=true) where {T} = +# ParamBox(Val(outSym), mapFunction, OneParam(inVar, inSym), genIndex(index), fill(canDiff)) + +# ParamBox(inVar::UMOTuple{NPMO, T}, outSym::Symbol, mapFunction::Function, +# inSym::NMOTuple{NPMO, Symbol}=ntuple(_->Symbol(IVsymSuffix, outSym), Val(NPMO)); +# index::Union{Int, Nothing}=nothing, canDiff::Bool=true) where {NPMO, T} = +# ParamBox(Val(outSym), mapFunction, +# OneParam.(inVar, inSym), genIndex(index), fill(canDiff)) + +const VPB{T} = Union{T, Array0D{T}, ParamBox{T}} -ParamBox(::Val{V}, pb::ParamBox{T}, index=genIndex(pb.canDiff[] ? pb.index[] : nothing), - canDiff::Array{Bool, 0}=copy(pb.canDiff)) where {V, T} = -ParamBox{T, V}(pb.map, pb.data[], index, canDiff) -ParamBox(inVar::Union{T, Array{T, 0}}, outSym::Symbol=:undef, - inSym::Symbol=Symbol(IVsymSuffix, outSym); - index::Union{Int, Nothing}=nothing, canDiff::Bool=false) where {T} = -ParamBox(Val(outSym), itself, - fillObj(inVar)=>inSym, genIndex(index), fill(canDiff)) +""" -ParamBox(inVar::Union{T, Array{T, 0}}, outSym::Symbol, mapFunction::Function, - inSym::Symbol=Symbol(IVsymSuffix, outSym); - index::Union{Int, Nothing}=nothing, canDiff::Bool=true) where {T} = -ParamBox(Val(outSym), mapFunction, - fillObj(inVar)=>inSym, genIndex(index), fill(canDiff)) + dataOf(pb::ParamBox{T}) where {T} -> Pair{Array{T, 0}, Symbol} -const VPB{T} = Union{T, Array{T, 0}, ParamBox{T}} +Return the `Pair` of the input variable and its symbol. +""" +@inline dataOf(pb::MonoParamBox) = pb.data +@inline dataOf(pb::PolyParamBox) = vcat(pb.data1, pb.data2, pb.data3) """ @@ -139,7 +261,7 @@ const VPB{T} = Union{T, Array{T, 0}, ParamBox{T}} Return the value of the input variable of `pb`. Equivalent to `pb[]`. """ -@inline inValOf(pb::ParamBox) = pb.data[][begin][] +@inline inValOf(pb::ParamBox) = valOf.(dataOf(pb)) """ @@ -148,7 +270,7 @@ Return the value of the input variable of `pb`. Equivalent to `pb[]`. Return the value of the output variable of `pb`. Equivalent to `pb()`. """ -@inline outValOf(pb::ParamBox) = (pb.map∘inValOf)(pb) +@inline outValOf(pb::ParamBox) = pb.map(inValOf(pb)) (pb::ParamBox)() = outValOf(pb) # (pb::ParamBox)() = Base.invokelatest(pb.map, pb.data[][begin][])::Float64 @@ -160,7 +282,7 @@ Return the value of the output variable of `pb`. Equivalent to `pb()`. Return the symbol of the input variable of `pb`. """ -@inline inSymOf(pb::ParamBox) = pb.data[][end] +@inline inSymOf(pb::ParamBox) = symOf.(dataOf(pb)) """ @@ -172,61 +294,57 @@ Return the symbol of the output variable of `pb`. @inline outSymOf(::ParamBox{<:Any, V}) where {V} = V -""" +# """ - isInSymEqual(pb::ParamBox, sym::Symbol) -> Bool +# isInSymEqual(pb::ParamBox, sym::Symbol) -> Bool -Return the Boolean value of whether the symbol of `pb`'s input variable equals `sym`. -""" -isInSymEqual(pb::ParamBox, sym::Symbol) = (dataOf(pb)[end] == sym) +# Return the Boolean value of whether the symbol of `pb`'s input variable equals `sym`. +# """ +# isInSymEqual(pb::ParamBox{<:Any, <:Any, <:Any, 1}, sym::TotTuple{Symbol}) = +# inSymOf(pb) == unzipObj(sym) +# isInSymEqual(pb::ParamBox{<:Any, <:Any, <:Any, NP}, sym::NTuple{NP, Symbol}) where {NP} = +# all(i==j for (i,j) in zip(inSymOfCore(pb), sym)) +# isInSymEqual(pb::ParamBox, sym::Symbol, i::Int) = symOf(pb.data[][i]) == sym -""" - isOutSymEqual(::ParamBox, sym::Symbol) -> Bool +# """ -Return the Boolean value of whether the symbol of `pb`'s output variable equals `sym`. -""" -isOutSymEqual(::ParamBox{<:Any, V}, sym::Symbol) where {V} = (V == sym) +# isOutSymEqual(::ParamBox, sym::Symbol) -> Bool + +# Return the Boolean value of whether the symbol of `pb`'s output variable equals `sym`. +# """ +# isOutSymEqual(::ParamBox{<:Any, V}, sym::Symbol) where {V} = (V == sym) """ - indVarOf(pb::ParamBox{T}) -> Pair{} + indParOf(pb::ParamBox{<:Any, <:Any, <:Any, NP}) -> Tuple{Vararg{SemiMutableParameter{T}, NP}} Return (the name and the value of) the independent variable tied to `pb`. Specifically, return the input variable stored in `pb` when `pb` is marked as differentiable; return the output variable of `pb` when `pb` is marked as non-differentiable. Thus, it is the variable `pb` represents to differentiate any (differentiable) function of [`ParamBox`](@ref)es. """ -function indVarOf(pb::ParamBox) - idx = numToSubs(pb.index[]) - if isDiffParam(pb) - Symbol(inSymOf(pb), idx) => inValOf(pb) - else - Symbol(outSymOf(pb), idx) => outValOf(pb) - end -end +indParOf(pb::ParamBox) = (unzipObj∘ifelse)(isDiffParam(pb), dataOf(pb), pb) +""" -getTypeParams(::Type{ParamBox{T, V, F}}) where {T, V, F} = (T, V, F) + indSymOf -getTypeParams(::T) where {T<:ParamBox} = getTypeParams(T) +""" +indSymOf(pb::ParamBox) = ifelse(isDiffParam(pb), outSymOf, inSymOf)(pb) -getFLevel(::Type{<:ParamBox{<:Any, <:Any, F}}) where {F} = getFLevel(F) - -getFLevel(::T) where {T<:ParamBox} = getFLevel(T) +getTypeParams(::Type{<:ParamBox{T, V, F}}) where {T, V, F} = (T, V, F) +getTypeParams(::T) where {T<:ParamBox} = getTypeParams(T) -""" - dataOf(pb::ParamBox{T}) where {T} -> Pair{Array{T, 0}, Symbol} +getFLevel(::Type{<:ParamBox{<:Any, <:Any, F}}) where {F} = getFLevel(F) -Return the `Pair` of the input variable and its symbol. -""" -@inline dataOf(pb::ParamBox) = pb.data[] +getFLevel(::T) where {T<:ParamBox} = getFLevel(T) """ @@ -246,8 +364,7 @@ Return a new `ParamBox` of which the input variable's value is equal to the outp variable's value of `pb`. """ outValCopy(pb::ParamBox{<:Any, V}) where {V} = -ParamBox(Val(V), itself, fill(pb())=>Symbol(IVsymSuffix, V)) - +genParamBox(Val(V), itself, OneParam(pb(), Symbol(IVsymSuffix, V))) """ @@ -255,7 +372,8 @@ ParamBox(Val(V), itself, fill(pb())=>Symbol(IVsymSuffix, V)) A shallow copy of the input `ParamBox`. """ -fullVarCopy(pb::ParamBox{<:Any, V}) where {V} = ParamBox(Val(V), pb, pb.index, pb.canDiff) +fullVarCopy(pb::ParamBox{<:Any, V}) where {V} = +genParamBox(Val(V), pb, pb.index[], pb.canDiff[]) """ @@ -314,7 +432,7 @@ end """ - changeMapping(pb::ParamBox{T, V}, mapFunction::Function=itself, outSym::Symbol=V; + changeMapping(pb::ParamBox{T, V}, mapFunction::Function, outSym::Symbol=V; canDiff::Union{Bool, Array{Bool, 0}}=isDiffParam(pb)) where {T, V} -> ParamBox{T, outSym} @@ -322,22 +440,35 @@ Change the mapping function of `pb`. The symbol of the output variable of the re `ParamBox` can be specified by `outSym`, and its differentiability is determined by `canDiff`. """ -function changeMapping(pb::ParamBox{T, V}, - mapFunction::Function=itself, outSym::Symbol=V; - canDiff::Union{Bool, Array{Bool, 0}}=isDiffParam(pb)) where {T, V} - canDiff = fillObj(canDiff) - ParamBox(Val(outSym), mapFunction, pb.data[], - genIndex( ifelse(canDiff[]==isDiffParam(pb)==true, pb.index[], nothing) ), - canDiff) +changeMapping(pb::ParamBox{T, V}, mapFunction::Function, outSym::Symbol=V; + canDiff::Union0D{Bool}=isDiffParam(pb)) where {T, V} = +genParamBox(Val(outSym), mapFunction, dataOf(pb), genIndex(), canDiff) + + +compareOneParam(op1::OneParam, op2::OneParam) = op1.data === op2.data + +compareParamBoxCore1(pb1::MonoParamBox{T}, pb2::MonoParamBox{T}) where {T} = +compareOneParam(dataOf(pb1)[], dataOf(pb2)[]) + +function compareParamBoxCore1(pb1::PloyParamBox{T}, pb2::PloyParamBox{T}) where {T} + dataPb1 = dataOf(pb1) + dataPb2 = dataOf(pb2) + if length(dataPb1) == length(dataPb2) + all( compareOneParam(op1, op2) for (op1, op2) in zip(dataPb1, dataPb2) ) + else + false + end end +# compareParamBoxCore1(pb1::ParamBox{<:Any, <:Any, <:Any, NP}, +# pb2::ParamBox{<:Any, <:Any, <:Any, NP}) where {NP} = +# all( compareOneParam(op1, op2) for (op1, op2) in zip(pb1.data[], pb2.data)) -compareParamBoxCore1(pb1::ParamBox, pb2::ParamBox) = -(pb1.data[][begin] === pb2.data[][begin]) +compareParamBoxCore1(pb1::ParamBox, pb2::ParamBox) = false -function compareParamBoxCore2(pb1::ParamBox{<:Any, V1, F1}, - pb2::ParamBox{<:Any, V2, F2}) where {V1, V2, F1, F2} - bl = V1==V2 && compareParamBoxCore1(pb1, pb2) +function compareParamBoxCore2(pb1::ParamBox{T1, V1, F1}, + pb2::ParamBox{T2, V2, F2}) where {T1, T2, V1, V2, F1, F2} + bl = T1==T2 && V1==V2 && compareParamBoxCore1(pb1, pb2) if FLevel(F1) == FLevel(F2) == IL bl else @@ -347,47 +478,41 @@ end @inline function compareParamBox(pb1::ParamBox, pb2::ParamBox) ifelse(( (bl=isDiffParam(pb1)) == isDiffParam(pb2) ), - ifelse( bl, - compareParamBoxCore1(pb1, pb2), - - compareParamBoxCore2(pb1, pb2) - ), - - false + (bl ? compareParamBoxCore1(pb1, pb2) : compareParamBoxCore2(pb1, pb2)), false ) end -function addParamBox(pb1::ParamBox{T, V, FL1}, pb2::ParamBox{T, V, FL2}, +function addParamBox(pb1::MonoParamBox{T, V, FL1}, pb2::MonoParamBox{T, V, FL2}, roundAtol::Real=nearestHalfOf(getAtolVal(T))) where {T, V, FL1, FL2} if isDiffParam(pb1) && compareParamBox(pb1, pb2) - ParamBox(Val(V), combinePF(+, pb1.map, pb2.map), - pb1.data[][begin]=>min(pb1.data[][end], pb2.data[][end]), - genIndex(nothing), fill(true)) + genParamBox(Val(V), combinePF(+, pb1.map, pb2.map), + pb1.data[][begin]=>min(pb1.data[][end], pb2.data[][end]), + genIndex(), fill(true)) else - ParamBox(Val(V), itself, + genParamBox(Val(V), itself, fill(roundToMultiOfStep(pb1() + pb2(), roundAtol))=>Symbol(IVsymSuffix, V)) end end -function mulParamBox(c::Number, pb::ParamBox{T, V}, +function mulParamBox(c::Number, pb::MonoParamBox{T, V}, roundAtol::Real=nearestHalfOf(getAtolVal(T))) where {T, V} if isDiffParam(pb) mapFunc = PF(pb.map, *, convert(T, roundToMultiOfStep(c, roundAtol))) - ParamBox(Val(V), mapFunc, pb.data[], genIndex(nothing), fill(true)) + ParamBox(Val(V), mapFunc, pb.data[], genIndex(), true) else ParamBox(Val(V), itself, fill(roundToMultiOfStep(T(pb()*c), roundAtol))=>Symbol(IVsymSuffix, V)) end end -function mulParamBox(pb1::ParamBox{T, V, FL1}, pb2::ParamBox{T, V, FL2}, +function mulParamBox(pb1::MonoParamBox{T, V, FL1}, pb2::MonoParamBox{T, V, FL2}, roundAtol::Real=nearestHalfOf(getAtolVal(T))) where {T, V, FL1, FL2} if isDiffParam(pb1) && compareParamBox(pb1, pb2) ParamBox(Val(V), combinePF(*, pb1.map, pb2.map), pb1.data[][begin]=>min(pb1.data[][end], pb2.data[][end]), - genIndex(nothing), fill(true)) + genIndex(), fill(true)) else ParamBox(Val(V), itself, fill(roundToMultiOfStep(pb1() * pb2(), roundAtol))=>Symbol(IVsymSuffix, V)) diff --git a/src/Quiqbox.jl b/src/Quiqbox.jl index 60aff3f0..f593850c 100644 --- a/src/Quiqbox.jl +++ b/src/Quiqbox.jl @@ -2,6 +2,7 @@ module Quiqbox include("AbstractTypes.jl") +include("Exception.jl") include("Tools.jl") include("FileIO.jl") diff --git a/src/Tools.jl b/src/Tools.jl index 44605b69..8f65c657 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -775,7 +775,8 @@ function genIndex(index::Int) genIndexCore(index) end -genIndex(index::Nothing) = genIndexCore(index) +genIndex() = genIndexCore(nothing) +genIndex(::Nothing) = genIndex() function genIndexCore(index) res = reshape(Union{Int, Nothing}[0], ()) |> collect @@ -797,14 +798,27 @@ function genNamedTupleC(name::Symbol, defaultVars::AbstractVector) end -fillObj(obj::Any) = fill(obj) +fillObj(@nospecialize(obj::Any)) = fill(obj) -fillObj(obj::Array{<:Any, 0}) = itself(obj) +fillObj(@nospecialize(obj::Array{<:Any, 0})) = itself(obj) -arrayToTuple(arr::AbstractArray) = Tuple(arr) +tupleObj(@nospecialize(obj::Any)) = (obj,) -arrayToTuple(tpl::Tuple) = itself(tpl) +tupleObj(@nospecialize(obj::Tuple{Any})) = itself(obj) + + +@inline unzipObj(obj::Tuple{T}) where {T} = obj[begin] + +unzipObj(@nospecialize(obj::Union{Tuple, Array})) = +ifelse(length(obj)==1, obj[begin], obj) + +unzipObj(@nospecialize(obj::Any)) = itself(obj) + + +arrayToTuple(@nospecialize(arr::AbstractArray)) = Tuple(arr) + +arrayToTuple(@nospecialize(tpl::Tuple)) = itself(tpl) genTupleCoords(::Type{T1}, diff --git a/test/AngularMomentum-test.jl b/test/AngularMomentum-test.jl new file mode 100644 index 00000000..985874e5 --- /dev/null +++ b/test/AngularMomentum-test.jl @@ -0,0 +1,62 @@ +using Quiqbox: numEps, hasApprox +using WignerSymbols +using Test + + +function testCG(j1, m1, j2, m2, j3) + + t = 2numEps(Float64) + + #= (m1-1, m2) (m1 , m2) =# + #= cg1___________cg3 =# + #= |\ | =# + #= | \ | =# + #= | \ | =# + #= | \ | =# + #= | \ | =# + #= |__________\| =# + #= cg4 cg2 =# + #= (m1-1,m2-1) (m1 ,m2-1) =# + + cg1_1 = CGcoeff(j1, m1-1, j2, m2, j3) + cg2_1 = CGcoeff(j1, m1, j2, m2-1, j3) + cg3_1 = genCGcoeff(+, cg1_1, cg2_1) + cg3_2 = genCGcoeff(+, cg2_1, cg1_1) + cg4_1 = genCGcoeff(-, cg1_1, cg2_1) + cg4_2 = genCGcoeff(-, cg2_1, cg1_1) + + @test cg3_1 == cg3_2 + @test cg3_1.m1 == m1 + @test cg3_1.m2 == m2 + @test cg4_1 == cg4_2 + @test cg4_1.m1 == m1-1 + @test cg4_1.m2 == m2-1 + @test typeof(cg3_1) == typeof(cg4_1) == CGcoeff{Float64, Int(2j1), Int(2j2), Int(2j3)} + @test hasApprox(cg3_1.coeff, getCGcoeff(cg3_1), getCGcoeff(j1, m1, j2, m2, j3), + clebschgordan(Float64, j1, m1, j2, m2, j3), atol=t) + @test hasApprox(cg4_1.coeff, getCGcoeff(cg4_1), getCGcoeff(j1, m1-1, j2, m2-1, j3), + clebschgordan(Float64, j1, m1-1, j2, m2-1, j3), atol=t) + + cg1_2 = genCGcoeff(-, cg3_1, cg2_1) + cg1_3 = genCGcoeff(-, cg2_1, cg3_1) + cg2_2 = genCGcoeff(-, cg3_1, cg1_1) + cg2_3 = genCGcoeff(-, cg1_1, cg3_1) + cg1_4 = genCGcoeff(+, cg4_1, cg2_1) + cg1_5 = genCGcoeff(+, cg2_1, cg4_1) + cg2_4 = genCGcoeff(+, cg4_1, cg1_1) + cg2_5 = genCGcoeff(+, cg1_1, cg4_1) + + @test cg1_1.m1 == cg1_2.m1 == cg1_3.m1 == cg1_4.m1 == cg1_5.m1 + @test cg1_1.m2 == cg1_2.m2 == cg1_3.m2 == cg1_4.m2 == cg1_5.m2 + @test hasApprox(cg1_1.coeff, cg1_2.coeff, cg1_3.coeff, cg1_4.coeff, cg1_5.coeff, atol=t) + + @test cg2_1.m1 == cg2_2.m1 == cg2_3.m1 == cg2_4.m1 == cg2_5.m1 + @test cg2_1.m2 == cg2_2.m2 == cg2_3.m2 == cg2_4.m2 == cg2_5.m2 + @test hasApprox(cg2_1.coeff, cg2_2.coeff, cg2_3.coeff, cg2_4.coeff, cg2_5.coeff, atol=t) +end + +j1m1j2m2j3_1 = (2, 2, 2, -1, 3) +testCG(j1m1j2m2j3_1...) + +j1m1j2m2j3_2 = (1.5, 0.5, 2, -1, 3.5) +testCG(j1m1j2m2j3_2...) diff --git a/test/runtests.jl b/test/runtests.jl index 8e06f918..e09d763b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,7 +32,7 @@ const SharedTestFunctions = true include("unit-tests/Integrals/Core-test.jl") include("unit-tests/Integrals/OneBody-test.jl") include("unit-tests/Integrals/TwoBody-test.jl") - Sys.islinux() && include("unit-tests/Integrals/Libcint-compare-tests.jl") + # Sys.islinux() && include("unit-tests/Integrals/Libcint-compare-tests.jl") include("unit-tests/Differentiation-test.jl") end println("$(unit2_1) test finished in $t2_1 seconds.\n") diff --git a/test/unit-tests/Basis-test.jl b/test/unit-tests/Basis-test.jl index 596c174e..b1efaf73 100644 --- a/test/unit-tests/Basis-test.jl +++ b/test/unit-tests/Basis-test.jl @@ -904,7 +904,7 @@ bf_gv7 = genBasisFunc(sp_gv1, (gf_gv4, gf_gv5)) pars_gv = markParams!(bf_gv7, true) pars_gv_t = markParams!([pb_gv1, pb_gv2, pb_gv3, gf_gv5.con]) @test pars_gv == pars_gv_t -@test (first∘indVarOf).(pars_gv) == [:I₁, :I₂, :I₃, :d₁] +@test indSymOf.(pars_gv) == [:I₁, :I₂, :I₃, :d₁] # function hasNormFactor getNormFactor diff --git a/test/unit-tests/Parameters-test.jl b/test/unit-tests/Parameters-test.jl index 8c2020a0..6fa63e0d 100644 --- a/test/unit-tests/Parameters-test.jl +++ b/test/unit-tests/Parameters-test.jl @@ -21,7 +21,7 @@ toggleDiff!(pb1) @test pb1.canDiff[] == false @test enableDiff!(pb1) @test pb1.canDiff[] == true -pb2 = changeMapping(pb1) +pb2 = changeMapping(pb1, itself) @test dataOf(pb1) === dataOf(pb2) === pb1.data[] pb2_2 = outValCopy(pb1) @test dataOf(pb1) == dataOf(pb2_2) diff --git a/test/unit-tests/Tools-test.jl b/test/unit-tests/Tools-test.jl index 9032a25b..0827893b 100644 --- a/test/unit-tests/Tools-test.jl +++ b/test/unit-tests/Tools-test.jl @@ -177,7 +177,7 @@ d = (2,4,5,1,3) # function genIndex @test genIndex(1) == fill(1) -@test genIndex(nothing) == fill(nothing) +@test genIndex() == fill(nothing) # function fillObj From df37e2a53ecee9bc6aaad8aea8786e6789db259a Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Sat, 4 Nov 2023 16:23:46 -0700 Subject: [PATCH 02/22] Optimized performance. --- src/Integrals/Core.jl | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index aec45aa4..a9255b2c 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -10,7 +10,7 @@ using Base: OneTo, Iterators.product ## [DOI] 10.48550/arXiv.2007.12057 function genFγIntegrand(γ::Int, u::T) where {T} - function (x) + function (x::T) ( (x+1)/2 )^(2γ) * exp(-u * (x+1)^2 / 4) / 2 end end @@ -25,19 +25,18 @@ for ValI in ValInts[begin:end .<= 1000] end function F0(u::T) where {T} - ifelse(u < getAtolVal(T), - T(1), - begin - ur = sqrt(u) - T(πvals[0.5]) * erf(ur) / (2ur) - end - ) + if u < getAtolVal(T) + T(1) + else + ur = sqrt(u) + T(πvals[0.5]) * erf(ur) / (2ur) + end end function getGQN(u::T) where {T} u = abs(u) + getAtolVal(T) - res = getAtolDigits(T) + round(0.4u + 2inv(sqrt(u))) + 1 - (Int∘min)(res, typemax(Int) - 1) + res = getAtolDigits(T) + (Int∘round)(0.4u + 2inv(sqrt(u))) + 1 + min(res, typemax(Int) - 1) end function Fγ(γ::Int, u::T) where {T} From 1b6209392b05b4891de67aa23e6b75af745a82ae Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Fri, 29 Dec 2023 22:11:23 -0500 Subject: [PATCH 03/22] Reverted back to pre-param-box opt stage. --- Project.toml | 2 - src/AbstractTypes.jl | 46 +---- src/Basis.jl | 41 +++-- src/Integrals/Core.jl | 2 + src/Overload.jl | 27 +-- src/Parameters.jl | 325 +++++++++++----------------------- src/Tools.jl | 4 +- test/unit-tests/Basis-test.jl | 2 +- 8 files changed, 137 insertions(+), 312 deletions(-) diff --git a/Project.toml b/Project.toml index d0a3b537..c70a8e83 100644 --- a/Project.toml +++ b/Project.toml @@ -11,11 +11,9 @@ LBFGSB = "5be7bae1-8223-5378-bac3-9e7378a2f6e6" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" -ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" SPGBox = "bf97046b-3e66-4aa0-9aed-26efb7fac769" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" [compat] diff --git a/src/AbstractTypes.jl b/src/AbstractTypes.jl index df948e83..d2a6a2c3 100644 --- a/src/AbstractTypes.jl +++ b/src/AbstractTypes.jl @@ -39,10 +39,8 @@ abstract type ManyFermionDataBox{T, D} <: ImmutableDataBox{T} end abstract type HartreeFockintermediateData{T} <: MutableDataBox{T} end -abstract type VariableBox{T, ContainerT} <: SemiMutableParameter{T, ContainerT} end -abstract type ConfigBox{T, ContainerT, MethodT} <: MutableParameter{T, ContainerT} end - -abstract type ParamBox{T, V, F} <: VariableBox{T, ParamBox} end +abstract type DifferentiableParameter{DataT, ContainerT} <: SemiMutableParameter{DataT, ContainerT} end +abstract type ConfigBox{T, ContainerT, MethodT} <: MutableParameter{ContainerT, T} end abstract type HartreeFockFinalValue{T, HFT} <: AbstractHartreeFockFinalValue{T} end @@ -103,42 +101,4 @@ const NTupleOfSBN{BN, T, D} = NTuple{BN, SpatialBasis{T, D}} const StrOrSym = Union{String, Symbol} -const MatrixCol{T} = Vector{SubArray{T, 1, Matrix{T}, Tuple{Slice{OneTo{Int}}, Int}, true}} - -const Array0D{T} = Array{T, 0} - -const Union0D{T} = Union{T, Array0D{T}} - -const MTuple{M, T} = NTuple{M, Array0D{T}} # M: mutable - -const UTuple{U, T} = NTuple{U, Union0D{T}} # U: union - -const NMOTuple{NMO, T} = Tuple{T, Vararg{T, NMO}} - -const MMOTuple{MMO, T} = Tuple{Array0D{T}, Vararg{Array0D{T}, MMO}} - -const UMOTuple{UMO, T} = Tuple{Union0D{T}, Vararg{Union0D{T}, UMO}} - -const NMTTuple{NMO, T} = Tuple{T, T, Vararg{T, NMO}} - -const MMTTuple{MMO, T} = Tuple{Array0D{T}, Array0D{T}, Vararg{Array0D{T}, MMO}} - -const UMTTuple{UMO, T} = Tuple{Union0D{T}, Union0D{T}, Vararg{Union0D{T}, UMO}} - -const TorTuple{T} = Union{T, Tuple{T}} - -const TorTupleLong{T} = Union{T, Tuple{T, Vararg{T}}} - -const TorAbstractVec{T} = Union{T, AbstractVector{<:T}} - -const IntOrNone = Union{Int, Nothing} -const IntOrMiss = Union{Int, Missing} - -const MutableIndex = Array0D{IntOrNone} - - -const MonoParam{T} = OneParam{T} - -const DualParam{T} = Tuple{OneParam{T}, OneParam{T}} - -const MoreParam{T} = Tuple{OneParam{T}, OneParam{T}, OneParam{T}, Vector{OneParam{T}}} \ No newline at end of file +const MatrixCol{T} = Vector{SubArray{T, 1, Matrix{T}, Tuple{Slice{OneTo{Int}}, Int}, true}} \ No newline at end of file diff --git a/src/Basis.jl b/src/Basis.jl index 61fa1e44..8b78c6e3 100644 --- a/src/Basis.jl +++ b/src/Basis.jl @@ -68,25 +68,24 @@ getTypeParams(::GaussFunc{T, Fxpn, Fcon}) where {T, Fxpn, Fcon} = (T, Fxpn, Fcon {T<:AbstractFloat} -> ParamBox{T, :$(xpnSym)} -Construct an exponent coefficient given a value or variable (and its symbol without the -index). `mapFunction`, `canDiff`, and `inSym` work the same way as in a general constructor -of a [`ParamBox`](@ref). +Construct an exponent coefficient given a value or variable (with its symbol). +`mapFunction`, `canDiff`, and `inSym` work the same way as in a general constructor of a +[`ParamBox`](@ref). """ -genExponent(eData::TorTuple{OneParam{T}}, mapFunction::Function=itself; +genExponent(eData::Pair{Array{T, 0}, Symbol}, mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true)) where {T<:AbstractFloat} = -ParamBox(Val(xpnSym), mapFunction, eData, genIndex(), fill(canDiff)) +ParamBox(Val(xpnSym), mapFunction, eData, genIndex(nothing), fill(canDiff)) -genExponent(eData::TorTuple{Union0D{T}}, mapFunction::Function=itself; +genExponent(e::Array{T, 0}, mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true), - inSym::TorTuple{Symbol}=(xpnIVsym,)) where {T<:AbstractFloat} = -genExponent(OneParam.(tupleObj(eData), inSym), mapFunction; canDiff) + inSym::Symbol=xpnIVsym) where {T<:AbstractFloat} = +genExponent(e=>inSym, mapFunction; canDiff) -genExponent(eData::NMOTuple{NPMO, Union{Union0D{T}, OneParam{T}}}, - mapFunction::Function=itself; +genExponent(e::AbstractFloat, mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true), - inSym::Symbol=xpnIVsym) where {T<:AbstractFloat} = -genExponent(OneParam.(tupleObj(eData), inSym), mapFunction; canDiff) + inSym::Symbol=xpnIVsym) = +genExponent(fill(e)=>inSym, mapFunction; canDiff) """ @@ -115,15 +114,15 @@ genExponent(pb::ParamBox) = ParamBox(Val(xpnSym), pb) {T<:AbstractFloat} -> ParamBox{T, :$(conSym)} -Construct a contraction coefficient given a value or variable (and its symbol without the -index).`mapFunction`, `canDiff`, and `inSym` work the same way as in a general constructor -of a [`ParamBox`](@ref). +Construct a contraction coefficient given a value or variable (with its symbol). +`mapFunction`, `canDiff`, and `inSym` work the same way as in a general constructor of a +[`ParamBox`](@ref). """ genContraction(dData::Pair{Array{T, 0}, Symbol}, mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true)) where {T<:AbstractFloat} = -ParamBox(Val(conSym), mapFunction, dData, genIndex(), fill(canDiff)) +ParamBox(Val(conSym), mapFunction, dData, genIndex(nothing), fill(canDiff)) genContraction(d::Array{T, 0}, mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true), @@ -269,7 +268,7 @@ genSpatialPoint(compData::Pair{Array{T, 0}, Symbol}, compIndex::Int, mapFunction::Function=itself; canDiff::Bool=ifelse(FLevel(mapFunction)==IL, false, true)) where {T<:AbstractFloat} = -ParamBox(Val(SpatialParamSyms[compIndex]), mapFunction, compData, genIndex(), +ParamBox(Val(SpatialParamSyms[compIndex]), mapFunction, compData, genIndex(nothing), fill(canDiff)) genSpatialPoint(comp::Array{T, 0}, compIndex::Int, @@ -1869,9 +1868,9 @@ end Return the parameter(s) stored in the input container. If `symbol` is set to `missing`, then return all the parameter(s). If it's set to a `Symbol` tied to a parameter, for example `:α₁`, the function will match any `pb::`[`ParamBox`](@ref) such that -[`indSymOf`](@ref)`(pb) == :α₁`. If it's set to a `Symbol` without any subscript, +[`indVarOf`](@ref)`(pb)[begin] == :α₁`. If it's set to a `Symbol` without any subscript, for example `:α`, the function will match it with all the `pb`s such that -`string(indSymOf(pb))` contains `'α'`. If the first argument is a collection, its +`string(indVarOf(pb)[begin])` contains `'α'`. If the first argument is a collection, its entries must be `ParamBox` containers. """ getParams(pb::ParamBox, symbol::Union{Symbol, Missing}=missing) = @@ -1909,7 +1908,7 @@ getParams(@nospecialize(cs::Tuple{QuiqboxContainer, Vararg{QuiqboxContainer}}), getParams(lazyCollect(cs), symbol) paramFilter(pb::ParamBox, sym::Union{Symbol, Missing}) = -sym isa Missing || inSymbol(sym, indSymOf(pb)) +sym isa Missing || inSymbol(sym, indVarOf(pb)[begin]) copyParContainer(pb::ParamBox, @@ -1990,7 +1989,7 @@ independent variable (in other words, considered identical in the differentiatio procedure) will be marked with same index. `filterMapping` determines whether filtering out (i.e. not return) the extra `ParamBox`es that hold the same independent variables but represent different output variables. The independent variable held by a `ParamBox` can be -inspected using [`indSymOf`](@ref). +inspected using [`indVarOf`](@ref). """ markParams!(b::Union{AbstractVector{T}, T}, filterMapping::Bool=false) where {T<:ParameterizedContainer} = diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index a9255b2c..83017cd6 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -17,6 +17,8 @@ end @generated function FγCore(γ::Int, u::T, ::Val{GQN}) where {T, GQN} GQnodes, GQweights = gausslegendre(GQN) + GQnodes = convert(Vector{T}, GQnodes) + GQweights = convert(Vector{T}, GQweights) return :(dot($GQweights, genFγIntegrand(γ, u).($GQnodes))) end diff --git a/src/Overload.jl b/src/Overload.jl index dddb097c..4d7c8ac1 100644 --- a/src/Overload.jl +++ b/src/Overload.jl @@ -1,13 +1,9 @@ # Julia internal methods overload. import Base: == -==(op1::OneParam, op2::OneParam) = op1.data == op2.data && - op1.label == op2.label && - op1.index[] === op2.index[] - ==(pb1::ParamBox, pb2::ParamBox) = (pb1.data == pb2.data && isDiffParam(pb1) == isDiffParam(pb2) && - pb1.index[] === pb2.index[] && + pb1.index[] == pb2.index[] && pb1.map == pb2.map) ==(cp1::SpatialPoint, cp2::SpatialPoint) = (cp1.param == cp2.param) @@ -89,7 +85,7 @@ end function printIndVarCore(io::IO, pb::ParamBox) - printstyled(io, "$(indSymOf(pb))", color=:cyan) + printstyled(io, "$(indVarOf(pb)[begin])", color=:cyan) end function printIndVar(io::IO, @nospecialize(pb::ParamBox)) @@ -304,14 +300,12 @@ import Base: * ## Iteration Interface import Base: iterate, size, length, eltype - -iterate(pb::ParamBox) = iterate(pb.data[]) -iterate(pb::ParamBox, i) = iterate(pb.data[], i) -size(::ParamBox{<:Any, <:Any, <:Any, 1}) = () -size(::ParamBox{<:Any, <:Any, <:Any, NP}) where {NP} = (NP,) -size(pb::ParamBox, d::Integer) = size(pb.data[], d) -length(::ParamBox{<:Any, <:Any, <:Any, NP}) where {NP} = NP +iterate(pb::ParamBox) = (pb.data[][begin][], nothing) +iterate(@nospecialize(_::ParamBox), _) = nothing +size(::ParamBox) = () +length(::ParamBox) = 1 eltype(::ParamBox{T}) where {T} = T +size(pb::ParamBox, d::Integer) = size(pb.data[][begin][], d) iterate(sp::SpatialPoint) = iterate(sp.param) iterate(sp::SpatialPoint, state) = iterate(sp.param, state) @@ -356,13 +350,10 @@ end ## Indexing Interface import Base: getindex, setindex!, firstindex, lastindex, eachindex, axes -getindex(op::OneParam) = op.data[] -getindex(pb::ParamBox) = inValOf(pb) +getindex(pb::ParamBox) = pb.data[][begin][] getindex(pb::ParamBox, ::Val{:first}) = getindex(pb) getindex(pb::ParamBox, ::Val{:last}) = getindex(pb) -setindex!(op::OneParam, d) = begin op.data[] = d end -setindex!(pb::ParamBox{<:Any, <:Any, <:Any, 1}, d) = begin pb.data[][begin] = d end -setindex!(pb::ParamBox, d, index) = begin pb.data[][index] = d end +setindex!(pb::ParamBox, d) = begin pb.data[][begin][] = d end firstindex(@nospecialize(_::ParamBox)) = Val(:first) lastindex(@nospecialize(_::ParamBox)) = Val(:last) axes(@nospecialize(_::ParamBox)) = () diff --git a/src/Parameters.jl b/src/Parameters.jl index b813558e..32b03638 100644 --- a/src/Parameters.jl +++ b/src/Parameters.jl @@ -1,33 +1,6 @@ export ParamBox, inValOf, outValOf, inSymOf, outSymOf, isInSymEqual, isOutSymEqual, - indParOf, indSymOf, dataOf, mapOf, outValCopy, fullVarCopy, enableDiff!, - disableDiff!, isDiffParam, toggleDiff!, changeMapping - -using StructArrays - -struct OneParam{T} <: VariableBox{T, OneParam} - label::Symbol - data::Array0D{T} - index::MutableIndex -end - -OneParam(label::Symbol, data::Union0D{T}, index::IntOrNone=nothing) where {T} = -OneParam{T}(label, fillObj(data), genIndex(index)) - -OneParam(labelToData::Pair{Symbol, <:Union0D{T}}, index::IntOrNone=nothing) where {T} = -OneParam{T}(labelToData[begin], fillObj(labelToData[end]), genIndex(index)) - -OneParam(op::OneParam) = itself(op) - -const OPessential{T} = Union{OneParam{T}, Pair{Symbol, <:Union0D{T}}} - -const OneParamStructVec{T} = - StructVector{OneParam{T}, NamedTuple{(:label, :data, :index), - Tuple{Vector{Symbol}, Vector{Array0D{T}}, Vector{MutableIndex}}}, Int} - - -valOf(par::OneParam) = par.data[] -symOf(par::OneParam) = Symbol(par.label, numToSubs(par.index[])) -# indexOf(par::OneParam) = par.index[] + indVarOf, dataOf, mapOf, outValCopy, fullVarCopy, enableDiff!, disableDiff!, + isDiffParam, toggleDiff!, changeMapping """ @@ -75,8 +48,7 @@ stored. If the latter is the type of `data`, then it will directly used to const `outSym::Symbol`: The symbol of the output variable represented by the constructed `ParamBox`. It's equal to the type parameter `V` of the constructed `ParamBox`. -`inSym::Symbol`: The symbol of the input variable (without the index) held by the -constructed `ParamBox`. +`inSym::Symbol`: The symbol of the input variable held by the constructed `ParamBox`. `mapFunction::Function`: The mapping (`mapFunction(::T)->T`) of the input variable, which will be assigned to the field `.map`. When `mapFunction` is not provided, `.map` is set to @@ -116,143 +88,49 @@ parameter functional (e.g., the Hartree–Fock energy). However, the derivative to the stored input variable can also be computed to when the `ParamBox` is marked as differentiable. """ -struct MonoParamBox{T, V, F<:Function} <: ParamBox{T, V, F} - data::Array0D{OneParam{T}} - map::Union{F, DI{F}} - index::MutableIndex - canDiff::Array0D{Bool} - - function MonoParamBox{T, V}(mapFunction::F, data::OneParam{T}, - index::Array0D{<:IntOrNone}, canDiff::Array0D{Bool}) where - {T, V, F<:Function} - checkFuncReturn(mapFunction, :mapFunction, T, T) - new{T, V, dressOf(F)}(fill(data), mapFunction, index, canDiff) - end - - MonoParamBox{T, V}(::IF, data::OneParam{T}, index, canDiff) where {T, V} = - new{T, V, iT}(fill(data), itself, index, canDiff) -end - -const DefaultMPBoxMap = itself - -struct PolyParamBox{T, V, F<:Function} <: ParamBox{T, V, F} - data1::Array0D{OneParam{T}} - data2::Array0D{OneParam{T}} - data3::OneParamStructVec{T} +struct ParamBox{T, V, F<:Function} <: DifferentiableParameter{T, ParamBox} + data::Array{Pair{Array{T, 0}, Symbol}, 0} map::Union{F, DI{F}} - index::MutableIndex - canDiff::Array0D{Bool} - - function PolyParamBox{T, V}(mapFunction::F, data1::OneParam{T}, data2::OneParam{T}, - data3::AbstractVector{OneParam{T}}, - index::Array0D{<:IntOrNone}, canDiff::Array0D{Bool}) where - {T, V, F<:Function} - checkFuncReturn(mapFunction, :mapFunction, Vector{T}, T) - new{T, V, dressOf(F)}(fill(data1), fill(data2) StructArray(data3), - f, index, canDiff) + canDiff::Array{Bool, 0} + index::Array{Union{Int, Nothing}, 0} + + function ParamBox{T, V}(f::F, data::Pair{Array{T, 0}, Symbol}, index, canDiff) where + {T, V, F<:Function} + Base.return_types(f, (T,))[1] == T || + throw(AssertionError("The mapping function `f`: `$(f)` should return the same "* + "data type as its input argument.")) + new{T, V, dressOf(F)}(fill(data), f, canDiff, index) end -end -const DefaultPPBoxMap = sum - -genParamBox(::Val{V}, mapFunction::F, data::OneParam{T}, - index::Union0D{<:IntOrNone}=genIndex(), - canDiff::Union0D{Bool}=getFLevel(F)>0) where {T, V, F<:Function} = -MonoParamBox{T, V}(mapFunction, data, fillObj(index), fillObj(canDiff)) - -genParamBox(::Val{V}, mapFunction::F, data1::OneParam{T}, data2::OneParam{T}, - data3::AbstractVector{OneParam{T}}, - index::Union0D{<:IntOrNone}=genIndex(), - canDiff::Union0D{Bool}=false) where {T, V, F<:Function} = -PolyParamBox{T, V}(mapFunction, data1, data2, data3, fillObj(index), fillObj(canDiff)) - -function genParamBox(::Val{V}, mapFunction::F, data::AbstractVector{OneParam{T}}, - index::Union0D{<:IntOrNone}=genIndex(), - canDiff::Union0D{Bool}=getFLevel(F)>0) where {T, V, F<:Function} - isLenOne = checkCollectionMinLen(data, :data, 1) - if isLenOne - genParamBox(Val(V), mapFunction, data[], index, canDiff) - else - genParamBox(Val(V), mapFunction, data[begin], data[begin+1], data[begin+2:end], - index, canDiff) - end + ParamBox{T, V}(data::Pair{Array{T, 0}, Symbol}, index, canDiff) where {T, V} = + new{T, V, iT}(fill(data), itself, canDiff, index) end -genParamBox(::Val{V}, pb::MonoParamBox{T}, index::Union0D{<:IntOrNone}=genIndex(), - canDiff::Array0D{Bool}=copy(pb.canDiff)) where {V, T} = -genParamBox(Val(V), pb.map, pb.data[], index, canDiff) - -genParamBox(::Val{V}, pb::PolyParamBox{T}, index::Union0D{<:IntOrNone}=genIndex(), - canDiff::Array0D{Bool}=copy(pb.canDiff)) where {V, T} = -genParamBox(Val(V), pb.map, pb.data1[], pb.data2[], pb.data3, index, canDiff) - -genParamBox(data::OPessential{T}, outSym::Symbol=:undef, mapFunction::F=DefaultMPBoxMap; - index::IntOrNone=nothing, - canDiff::Bool=getFLevel(F)>0) where {T, F<:Function} = -genParamBox(Val(outSym), mapFunction, OneParam(data), index, fill(canDiff)) - -function genParamBox(data::AbstractVector{<:Union{Union0D{T}, OPessential{T}}}, - outSym::Symbol=:undef, - mapFunction::F=ifelse(checkCollectionMinLen(data, :data, 1), - DefaultMPBoxMap, DefaultPPBoxMap); - defaultDataLabel::Symbol=Symbol(IVsymSuffix, outSym), - index::IntOrNone=nothing, - canDiff::Bool=getFLevel(F)>0) where {T, F<:Function} - data = map(data) do item - item isa Union0D ? OneParam(defaultDataLabel, item) : OneParam(item) - end - genParamBox(Val(outSym), mapFunction, data, index, canDiff) -end +ParamBox(::Val{V}, mapFunction::F, data::Pair{Array{T, 0}, Symbol}, + index=genIndex(nothing), canDiff=fill(true)) where {V, F<:Function, T} = +ParamBox{T, V}(mapFunction, data, index, canDiff) -# function ParamBox(data::Union0D{T}, outSym::Symbol=:undef, mapFunction::F=itself; -# defaultDataLabel::Symbol=Symbol(IVsymSuffix, outSym), -# index::IntOrNone=nothing, canDiff::Bool=getFLevel(F)>0) where -# {T, F<:Function} -# data = map(data|>tupleObj) do item -# item isa Union0D ? OneParam(defaultDataLabel, item) : item -# end -# ParamBox(Val(outSym), mapFunction, OneParam.(tupleObj(data)), -# genIndex(index), fill(canDiff)) -# end - -# function ParamBox(data::TorTupleLong{Union{Union0D{T}, OPessential{T}}}, -# outSym::Symbol=:undef, mapFunction::Function=itself; -# defaultDataLabel::Symbol=Symbol(IVsymSuffix, outSym), -# index::Union{Int, Nothing}=nothing, canDiff::Bool=false) -# data = map(data|>tupleObj) do item -# item isa Union0D ? OneParam(defaultDataLabel, item) : item -# end -# ParamBox(Val(outSym), mapFunction, OneParam.(tupleObj(data)), -# genIndex(index), fill(canDiff)) -# end - -# ParamBox(inVar::Union0D{T}, outSym::Symbol=:undef, -# inSym::Symbol=(Symbol(IVsymSuffix, outSym),); -# index::Union{Int, Nothing}=nothing, canDiff::Bool=false) where {T} = -# ParamBox(Val(outSym), itself, OneParam(inVar, inSym), genIndex(index), fill(canDiff)) - -# ParamBox(inVar::Union0D{T}, outSym::Symbol, mapFunction::Function, -# inSym::Symbol=Symbol(IVsymSuffix, outSym); -# index::Union{Int, Nothing}=nothing, canDiff::Bool=true) where {T} = -# ParamBox(Val(outSym), mapFunction, OneParam(inVar, inSym), genIndex(index), fill(canDiff)) - -# ParamBox(inVar::UMOTuple{NPMO, T}, outSym::Symbol, mapFunction::Function, -# inSym::NMOTuple{NPMO, Symbol}=ntuple(_->Symbol(IVsymSuffix, outSym), Val(NPMO)); -# index::Union{Int, Nothing}=nothing, canDiff::Bool=true) where {NPMO, T} = -# ParamBox(Val(outSym), mapFunction, -# OneParam.(inVar, inSym), genIndex(index), fill(canDiff)) - -const VPB{T} = Union{T, Array0D{T}, ParamBox{T}} +ParamBox(::Val{V}, ::IF, data::Pair{Array{T, 0}, Symbol}, index=genIndex(nothing), + canDiff=fill(false)) where {V, T} = +ParamBox{T, V}(data, index, canDiff) +ParamBox(::Val{V}, pb::ParamBox{T}, index=genIndex(pb.canDiff[] ? pb.index[] : nothing), + canDiff::Array{Bool, 0}=copy(pb.canDiff)) where {V, T} = +ParamBox{T, V}(pb.map, pb.data[], index, canDiff) -""" +ParamBox(inVar::Union{T, Array{T, 0}}, outSym::Symbol=:undef, + inSym::Symbol=Symbol(IVsymSuffix, outSym); + index::Union{Int, Nothing}=nothing, canDiff::Bool=false) where {T} = +ParamBox(Val(outSym), itself, + fillObj(inVar)=>inSym, genIndex(index), fill(canDiff)) - dataOf(pb::ParamBox{T}) where {T} -> Pair{Array{T, 0}, Symbol} +ParamBox(inVar::Union{T, Array{T, 0}}, outSym::Symbol, mapFunction::Function, + inSym::Symbol=Symbol(IVsymSuffix, outSym); + index::Union{Int, Nothing}=nothing, canDiff::Bool=true) where {T} = +ParamBox(Val(outSym), mapFunction, + fillObj(inVar)=>inSym, genIndex(index), fill(canDiff)) -Return the `Pair` of the input variable and its symbol. -""" -@inline dataOf(pb::MonoParamBox) = pb.data -@inline dataOf(pb::PolyParamBox) = vcat(pb.data1, pb.data2, pb.data3) +const VPB{T} = Union{T, Array{T, 0}, ParamBox{T}} """ @@ -261,7 +139,7 @@ Return the `Pair` of the input variable and its symbol. Return the value of the input variable of `pb`. Equivalent to `pb[]`. """ -@inline inValOf(pb::ParamBox) = valOf.(dataOf(pb)) +@inline inValOf(pb::ParamBox) = pb.data[][begin][] """ @@ -270,7 +148,7 @@ Return the value of the input variable of `pb`. Equivalent to `pb[]`. Return the value of the output variable of `pb`. Equivalent to `pb()`. """ -@inline outValOf(pb::ParamBox) = pb.map(inValOf(pb)) +@inline outValOf(pb::ParamBox) = (pb.map∘inValOf)(pb) (pb::ParamBox)() = outValOf(pb) # (pb::ParamBox)() = Base.invokelatest(pb.map, pb.data[][begin][])::Float64 @@ -282,7 +160,7 @@ Return the value of the output variable of `pb`. Equivalent to `pb()`. Return the symbol of the input variable of `pb`. """ -@inline inSymOf(pb::ParamBox) = symOf.(dataOf(pb)) +@inline inSymOf(pb::ParamBox) = pb.data[][end] """ @@ -294,50 +172,45 @@ Return the symbol of the output variable of `pb`. @inline outSymOf(::ParamBox{<:Any, V}) where {V} = V -# """ - -# isInSymEqual(pb::ParamBox, sym::Symbol) -> Bool +""" -# Return the Boolean value of whether the symbol of `pb`'s input variable equals `sym`. -# """ -# isInSymEqual(pb::ParamBox{<:Any, <:Any, <:Any, 1}, sym::TotTuple{Symbol}) = -# inSymOf(pb) == unzipObj(sym) + isInSymEqual(pb::ParamBox, sym::Symbol) -> Bool -# isInSymEqual(pb::ParamBox{<:Any, <:Any, <:Any, NP}, sym::NTuple{NP, Symbol}) where {NP} = -# all(i==j for (i,j) in zip(inSymOfCore(pb), sym)) +Return the Boolean value of whether the symbol of `pb`'s input variable equals `sym`. +""" +isInSymEqual(pb::ParamBox, sym::Symbol) = (dataOf(pb)[end] == sym) -# isInSymEqual(pb::ParamBox, sym::Symbol, i::Int) = symOf(pb.data[][i]) == sym -# """ +""" -# isOutSymEqual(::ParamBox, sym::Symbol) -> Bool + isOutSymEqual(::ParamBox, sym::Symbol) -> Bool -# Return the Boolean value of whether the symbol of `pb`'s output variable equals `sym`. -# """ -# isOutSymEqual(::ParamBox{<:Any, V}, sym::Symbol) where {V} = (V == sym) +Return the Boolean value of whether the symbol of `pb`'s output variable equals `sym`. +""" +isOutSymEqual(::ParamBox{<:Any, V}, sym::Symbol) where {V} = (V == sym) """ - indParOf(pb::ParamBox{<:Any, <:Any, <:Any, NP}) -> Tuple{Vararg{SemiMutableParameter{T}, NP}} + indVarOf(pb::ParamBox{T}) -> Pair{} Return (the name and the value of) the independent variable tied to `pb`. Specifically, return the input variable stored in `pb` when `pb` is marked as differentiable; return the output variable of `pb` when `pb` is marked as non-differentiable. Thus, it is the variable `pb` represents to differentiate any (differentiable) function of [`ParamBox`](@ref)es. """ -indParOf(pb::ParamBox) = (unzipObj∘ifelse)(isDiffParam(pb), dataOf(pb), pb) - -""" - - indSymOf - -""" -indSymOf(pb::ParamBox) = ifelse(isDiffParam(pb), outSymOf, inSymOf)(pb) +function indVarOf(pb::ParamBox) + idx = numToSubs(pb.index[]) + if isDiffParam(pb) + Symbol(inSymOf(pb), idx) => inValOf(pb) + else + Symbol(outSymOf(pb), idx) => outValOf(pb) + end +end -getTypeParams(::Type{<:ParamBox{T, V, F}}) where {T, V, F} = (T, V, F) +getTypeParams(::Type{ParamBox{T, V, F}}) where {T, V, F} = (T, V, F) getTypeParams(::T) where {T<:ParamBox} = getTypeParams(T) @@ -347,6 +220,15 @@ getFLevel(::Type{<:ParamBox{<:Any, <:Any, F}}) where {F} = getFLevel(F) getFLevel(::T) where {T<:ParamBox} = getFLevel(T) +""" + + dataOf(pb::ParamBox{T}) where {T} -> Pair{Array{T, 0}, Symbol} + +Return the `Pair` of the input variable and its symbol. +""" +@inline dataOf(pb::ParamBox) = pb.data[] + + """ mapOf(pb::ParamBox) -> Function @@ -364,7 +246,8 @@ Return a new `ParamBox` of which the input variable's value is equal to the outp variable's value of `pb`. """ outValCopy(pb::ParamBox{<:Any, V}) where {V} = -genParamBox(Val(V), itself, OneParam(pb(), Symbol(IVsymSuffix, V))) +ParamBox(Val(V), itself, fill(pb())=>Symbol(IVsymSuffix, V)) + """ @@ -372,8 +255,7 @@ genParamBox(Val(V), itself, OneParam(pb(), Symbol(IVsymSuffix, V))) A shallow copy of the input `ParamBox`. """ -fullVarCopy(pb::ParamBox{<:Any, V}) where {V} = -genParamBox(Val(V), pb, pb.index[], pb.canDiff[]) +fullVarCopy(pb::ParamBox{<:Any, V}) where {V} = ParamBox(Val(V), pb, pb.index, pb.canDiff) """ @@ -432,7 +314,7 @@ end """ - changeMapping(pb::ParamBox{T, V}, mapFunction::Function, outSym::Symbol=V; + changeMapping(pb::ParamBox{T, V}, mapFunction::Function=itself, outSym::Symbol=V; canDiff::Union{Bool, Array{Bool, 0}}=isDiffParam(pb)) where {T, V} -> ParamBox{T, outSym} @@ -440,35 +322,22 @@ Change the mapping function of `pb`. The symbol of the output variable of the re `ParamBox` can be specified by `outSym`, and its differentiability is determined by `canDiff`. """ -changeMapping(pb::ParamBox{T, V}, mapFunction::Function, outSym::Symbol=V; - canDiff::Union0D{Bool}=isDiffParam(pb)) where {T, V} = -genParamBox(Val(outSym), mapFunction, dataOf(pb), genIndex(), canDiff) - - -compareOneParam(op1::OneParam, op2::OneParam) = op1.data === op2.data - -compareParamBoxCore1(pb1::MonoParamBox{T}, pb2::MonoParamBox{T}) where {T} = -compareOneParam(dataOf(pb1)[], dataOf(pb2)[]) - -function compareParamBoxCore1(pb1::PloyParamBox{T}, pb2::PloyParamBox{T}) where {T} - dataPb1 = dataOf(pb1) - dataPb2 = dataOf(pb2) - if length(dataPb1) == length(dataPb2) - all( compareOneParam(op1, op2) for (op1, op2) in zip(dataPb1, dataPb2) ) - else - false - end +function changeMapping(pb::ParamBox{T, V}, + mapFunction::Function=itself, outSym::Symbol=V; + canDiff::Union{Bool, Array{Bool, 0}}=isDiffParam(pb)) where {T, V} + canDiff = fillObj(canDiff) + ParamBox(Val(outSym), mapFunction, pb.data[], + genIndex( ifelse(canDiff[]==isDiffParam(pb)==true, pb.index[], nothing) ), + canDiff) end -# compareParamBoxCore1(pb1::ParamBox{<:Any, <:Any, <:Any, NP}, -# pb2::ParamBox{<:Any, <:Any, <:Any, NP}) where {NP} = -# all( compareOneParam(op1, op2) for (op1, op2) in zip(pb1.data[], pb2.data)) -compareParamBoxCore1(pb1::ParamBox, pb2::ParamBox) = false +compareParamBoxCore1(pb1::ParamBox, pb2::ParamBox) = +(pb1.data[][begin] === pb2.data[][begin]) -function compareParamBoxCore2(pb1::ParamBox{T1, V1, F1}, - pb2::ParamBox{T2, V2, F2}) where {T1, T2, V1, V2, F1, F2} - bl = T1==T2 && V1==V2 && compareParamBoxCore1(pb1, pb2) +function compareParamBoxCore2(pb1::ParamBox{<:Any, V1, F1}, + pb2::ParamBox{<:Any, V2, F2}) where {V1, V2, F1, F2} + bl = V1==V2 && compareParamBoxCore1(pb1, pb2) if FLevel(F1) == FLevel(F2) == IL bl else @@ -478,41 +347,47 @@ end @inline function compareParamBox(pb1::ParamBox, pb2::ParamBox) ifelse(( (bl=isDiffParam(pb1)) == isDiffParam(pb2) ), - (bl ? compareParamBoxCore1(pb1, pb2) : compareParamBoxCore2(pb1, pb2)), false + ifelse( bl, + compareParamBoxCore1(pb1, pb2), + + compareParamBoxCore2(pb1, pb2) + ), + + false ) end -function addParamBox(pb1::MonoParamBox{T, V, FL1}, pb2::MonoParamBox{T, V, FL2}, +function addParamBox(pb1::ParamBox{T, V, FL1}, pb2::ParamBox{T, V, FL2}, roundAtol::Real=nearestHalfOf(getAtolVal(T))) where {T, V, FL1, FL2} if isDiffParam(pb1) && compareParamBox(pb1, pb2) - genParamBox(Val(V), combinePF(+, pb1.map, pb2.map), - pb1.data[][begin]=>min(pb1.data[][end], pb2.data[][end]), - genIndex(), fill(true)) + ParamBox(Val(V), combinePF(+, pb1.map, pb2.map), + pb1.data[][begin]=>min(pb1.data[][end], pb2.data[][end]), + genIndex(nothing), fill(true)) else - genParamBox(Val(V), itself, + ParamBox(Val(V), itself, fill(roundToMultiOfStep(pb1() + pb2(), roundAtol))=>Symbol(IVsymSuffix, V)) end end -function mulParamBox(c::Number, pb::MonoParamBox{T, V}, +function mulParamBox(c::Number, pb::ParamBox{T, V}, roundAtol::Real=nearestHalfOf(getAtolVal(T))) where {T, V} if isDiffParam(pb) mapFunc = PF(pb.map, *, convert(T, roundToMultiOfStep(c, roundAtol))) - ParamBox(Val(V), mapFunc, pb.data[], genIndex(), true) + ParamBox(Val(V), mapFunc, pb.data[], genIndex(nothing), fill(true)) else ParamBox(Val(V), itself, fill(roundToMultiOfStep(T(pb()*c), roundAtol))=>Symbol(IVsymSuffix, V)) end end -function mulParamBox(pb1::MonoParamBox{T, V, FL1}, pb2::MonoParamBox{T, V, FL2}, +function mulParamBox(pb1::ParamBox{T, V, FL1}, pb2::ParamBox{T, V, FL2}, roundAtol::Real=nearestHalfOf(getAtolVal(T))) where {T, V, FL1, FL2} if isDiffParam(pb1) && compareParamBox(pb1, pb2) ParamBox(Val(V), combinePF(*, pb1.map, pb2.map), pb1.data[][begin]=>min(pb1.data[][end], pb2.data[][end]), - genIndex(), fill(true)) + genIndex(nothing), fill(true)) else ParamBox(Val(V), itself, fill(roundToMultiOfStep(pb1() * pb2(), roundAtol))=>Symbol(IVsymSuffix, V)) diff --git a/src/Tools.jl b/src/Tools.jl index 8f65c657..e0f8fec9 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -775,8 +775,8 @@ function genIndex(index::Int) genIndexCore(index) end -genIndex() = genIndexCore(nothing) -genIndex(::Nothing) = genIndex() +genIndex(::Nothing) = genIndexCore(nothing) +genIndex() = genIndex(nothing) function genIndexCore(index) res = reshape(Union{Int, Nothing}[0], ()) |> collect diff --git a/test/unit-tests/Basis-test.jl b/test/unit-tests/Basis-test.jl index b1efaf73..596c174e 100644 --- a/test/unit-tests/Basis-test.jl +++ b/test/unit-tests/Basis-test.jl @@ -904,7 +904,7 @@ bf_gv7 = genBasisFunc(sp_gv1, (gf_gv4, gf_gv5)) pars_gv = markParams!(bf_gv7, true) pars_gv_t = markParams!([pb_gv1, pb_gv2, pb_gv3, gf_gv5.con]) @test pars_gv == pars_gv_t -@test indSymOf.(pars_gv) == [:I₁, :I₂, :I₃, :d₁] +@test (first∘indVarOf).(pars_gv) == [:I₁, :I₂, :I₃, :d₁] # function hasNormFactor getNormFactor From 88ab566b1a77fe8add41e8c838e0104914a5dee2 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Sat, 30 Dec 2023 02:58:52 -0500 Subject: [PATCH 04/22] Improved int func's performance. --- src/Integrals/Core.jl | 23 +++++++++++++++++++++-- 1 file changed, 21 insertions(+), 2 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 83017cd6..8de16b55 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -9,10 +9,29 @@ using Base: OneTo, Iterators.product ## [DOI] 10.1088/0143-0807/31/1/004 ## [DOI] 10.48550/arXiv.2007.12057 +# function genFγIntegrand(γ::Int, u::T) where {T} +# function (x::T) +# ( (x+1)/2 )^(2γ) * exp(-u * (x+1)^2 / 4) / 2 +# end +# end + +## Ineffective solution +# struct FγIntegrand{T} +# γ::Int +# u::T +# end +# (FγI::FγIntegrand{T})(x::T) where {T} = ( (x+1)/2 )^(2FγI.γ) * exp(-FγI.u * (x+1)^2 / 4) / 2 +# function genFγIntegrand(γ::Int, u::T) where {T} +# FγIntegrand(γ, u) +# end + function genFγIntegrand(γ::Int, u::T) where {T} - function (x::T) - ( (x+1)/2 )^(2γ) * exp(-u * (x+1)^2 / 4) / 2 + f = let γLoc=γ, uLoc=u + @inline function (x::T) # @inline does improve performance as of Julia 1.10.0 + ( (x+1)/2 )^(2γLoc) * exp(-uLoc * (x+1)^2 / 4) / 2 + end end + f end @generated function FγCore(γ::Int, u::T, ::Val{GQN}) where {T, GQN} From bb20aa62e1d2d030672a91fec276084d71942121 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Sat, 30 Dec 2023 02:59:51 -0500 Subject: [PATCH 05/22] Organized code. --- src/Integrals/Core.jl | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 8de16b55..4851e041 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -9,22 +9,6 @@ using Base: OneTo, Iterators.product ## [DOI] 10.1088/0143-0807/31/1/004 ## [DOI] 10.48550/arXiv.2007.12057 -# function genFγIntegrand(γ::Int, u::T) where {T} -# function (x::T) -# ( (x+1)/2 )^(2γ) * exp(-u * (x+1)^2 / 4) / 2 -# end -# end - -## Ineffective solution -# struct FγIntegrand{T} -# γ::Int -# u::T -# end -# (FγI::FγIntegrand{T})(x::T) where {T} = ( (x+1)/2 )^(2FγI.γ) * exp(-FγI.u * (x+1)^2 / 4) / 2 -# function genFγIntegrand(γ::Int, u::T) where {T} -# FγIntegrand(γ, u) -# end - function genFγIntegrand(γ::Int, u::T) where {T} f = let γLoc=γ, uLoc=u @inline function (x::T) # @inline does improve performance as of Julia 1.10.0 From e933ceb28acab9be2d67d41216d6f983e112a838 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Sun, 31 Dec 2023 02:55:17 -0500 Subject: [PATCH 06/22] Replaced generated func with cached func for type stability. --- Project.toml | 2 ++ src/Integrals/Core.jl | 53 ++++++++++++++++++++++++++++++++----------- src/Library.jl | 7 ------ test/runtests.jl | 2 +- 4 files changed, 43 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index c70a8e83..badac9ba 100644 --- a/Project.toml +++ b/Project.toml @@ -8,6 +8,7 @@ DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LBFGSB = "5be7bae1-8223-5378-bac3-9e7378a2f6e6" +LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -22,6 +23,7 @@ DoubleFloats = "^1.2" FastGaussQuadrature = "1" ForwardDiff = "^0.10.25" LBFGSB = "0.4" +LRUCache = "^1.6" LineSearches = "^7.1.1" LinearAlgebra = "1" Printf = "1" diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 4851e041..3d1c48c4 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -4,6 +4,7 @@ using SpecialFunctions: erf using FastGaussQuadrature: gausslegendre using LinearAlgebra: dot using Base: OneTo, Iterators.product +using LRUCache # Reference(s): ## [DOI] 10.1088/0143-0807/31/1/004 @@ -12,22 +13,48 @@ using Base: OneTo, Iterators.product function genFγIntegrand(γ::Int, u::T) where {T} f = let γLoc=γ, uLoc=u @inline function (x::T) # @inline does improve performance as of Julia 1.10.0 - ( (x+1)/2 )^(2γLoc) * exp(-uLoc * (x+1)^2 / 4) / 2 + ((x + 1) / 2)^(2γLoc) * exp(-uLoc * (x+1)^2 / 4) / 2 end end f end -@generated function FγCore(γ::Int, u::T, ::Val{GQN}) where {T, GQN} - GQnodes, GQweights = gausslegendre(GQN) +const NodesAndWeightsOfGQ = LRU{Int, Tuple{Vector{Float64}, Vector{Float64}}}(maxsize=200) + +const MaxGQpointNum = 10000 + +function FγCore(γ::Int, u::T, GQpointNum::Int) where {T} + GQpointNum = min(GQpointNum, MaxGQpointNum) + GQnodes, GQweights = get!(NodesAndWeightsOfGQ, GQpointNum) do + gausslegendre(GQpointNum) + end GQnodes = convert(Vector{T}, GQnodes) GQweights = convert(Vector{T}, GQweights) - return :(dot($GQweights, genFγIntegrand(γ, u).($GQnodes))) -end - -for ValI in ValInts[begin:end .<= 1000] - precompile(FγCore, (Int, Float64, ValI)) -end + dot(GQweights, genFγIntegrand(γ, u).(GQnodes)) +end + +## Multi-threaded version: too much overhead! +# const ThreadNum = Threads.nthreads() +# chopIter(len, startIdx=1) = Iterators.partition( startIdx:(startIdx+len-1), +# (1 + len ÷ 2ThreadNum) ) +# function FγCore(γ::Int, u::T, GQpointNum::Int) where {T} +# GQpointNum = min(GQpointNum, MaxGQpointNum) +# GQnodes, GQweights = get!(NodesAndWeightsOfGQ, GQpointNum) do +# gausslegendre(GQpointNum) +# end +# GQnodes = convert(Vector{T}, GQnodes) +# GQweights = convert(Vector{T}, GQweights) +# let f = genFγIntegrand(γ, u) +# tasks = map(chopIter(GQpointNum, firstindex(GQnodes))) do idx +# Threads.@spawn begin +# mapreduce(+, view(GQweights, idx), view(GQnodes, idx)) do weight, node +# weight * f(node) +# end +# end +# end +# mapreduce(fetch, +, tasks) +# end +# end function F0(u::T) where {T} if u < getAtolVal(T) @@ -38,17 +65,17 @@ function F0(u::T) where {T} end end -function getGQN(u::T) where {T} +function getGQpointNum(u::T) where {T} u = abs(u) + getAtolVal(T) - res = getAtolDigits(T) + (Int∘round)(0.4u + 2inv(sqrt(u))) + 1 + res = getAtolDigits(T) + (Int∘round)(0.4u + 2(inv∘sqrt)(u)) + 1 min(res, typemax(Int) - 1) end function Fγ(γ::Int, u::T) where {T} if u < getAtolVal(T) - (inv∘T∘muladd)(2, γ, 1) + (T∘inv∘muladd)(2, γ, 1) else - FγCore(γ, u, (getValI∘getGQN)(u)) + FγCore(γ, u, getGQpointNum(u)) end end diff --git a/src/Library.jl b/src/Library.jl index 02e5b7f6..c9e2347f 100644 --- a/src/Library.jl +++ b/src/Library.jl @@ -247,13 +247,6 @@ function checkBSList(;printInfo::Bool=false) end -const ValInts = Val.(collect(0:5000)) - -function getValI(i::Int) - i = abs(i) - i < length(ValInts) ? ValInts[i+1] : ValInts[end] -end - const πvals = Dict([-0.75, 0.5, 2.5] .=> big(π).^[-0.75, 0.5, 2.5]) const DefaultDigits = 10 diff --git a/test/runtests.jl b/test/runtests.jl index e09d763b..8e06f918 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,7 +32,7 @@ const SharedTestFunctions = true include("unit-tests/Integrals/Core-test.jl") include("unit-tests/Integrals/OneBody-test.jl") include("unit-tests/Integrals/TwoBody-test.jl") - # Sys.islinux() && include("unit-tests/Integrals/Libcint-compare-tests.jl") + Sys.islinux() && include("unit-tests/Integrals/Libcint-compare-tests.jl") include("unit-tests/Differentiation-test.jl") end println("$(unit2_1) test finished in $t2_1 seconds.\n") From 9d80c5b95346b661ff55b50440079d54a8099095 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 1 Jan 2024 11:16:53 -0500 Subject: [PATCH 07/22] Applied `LazyArrays` to reduce memory allocations. --- Project.toml | 2 ++ src/Integrals/Core.jl | 45 +++++++++++++------------------------------ src/Tools.jl | 20 +++++++++++++++++++ 3 files changed, 35 insertions(+), 32 deletions(-) diff --git a/Project.toml b/Project.toml index badac9ba..02fe0e25 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LBFGSB = "5be7bae1-8223-5378-bac3-9e7378a2f6e6" LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" +LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -24,6 +25,7 @@ FastGaussQuadrature = "1" ForwardDiff = "^0.10.25" LBFGSB = "0.4" LRUCache = "^1.6" +LazyArrays = "^1.8" LineSearches = "^7.1.1" LinearAlgebra = "1" Printf = "1" diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 3d1c48c4..9385fdcb 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -5,6 +5,7 @@ using FastGaussQuadrature: gausslegendre using LinearAlgebra: dot using Base: OneTo, Iterators.product using LRUCache +using LazyArrays: BroadcastArray # Reference(s): ## [DOI] 10.1088/0143-0807/31/1/004 @@ -26,35 +27,13 @@ const MaxGQpointNum = 10000 function FγCore(γ::Int, u::T, GQpointNum::Int) where {T} GQpointNum = min(GQpointNum, MaxGQpointNum) GQnodes, GQweights = get!(NodesAndWeightsOfGQ, GQpointNum) do - gausslegendre(GQpointNum) + gausslegendre(GQpointNum) end GQnodes = convert(Vector{T}, GQnodes) GQweights = convert(Vector{T}, GQweights) - dot(GQweights, genFγIntegrand(γ, u).(GQnodes)) -end - -## Multi-threaded version: too much overhead! -# const ThreadNum = Threads.nthreads() -# chopIter(len, startIdx=1) = Iterators.partition( startIdx:(startIdx+len-1), -# (1 + len ÷ 2ThreadNum) ) -# function FγCore(γ::Int, u::T, GQpointNum::Int) where {T} -# GQpointNum = min(GQpointNum, MaxGQpointNum) -# GQnodes, GQweights = get!(NodesAndWeightsOfGQ, GQpointNum) do -# gausslegendre(GQpointNum) -# end -# GQnodes = convert(Vector{T}, GQnodes) -# GQweights = convert(Vector{T}, GQweights) -# let f = genFγIntegrand(γ, u) -# tasks = map(chopIter(GQpointNum, firstindex(GQnodes))) do idx -# Threads.@spawn begin -# mapreduce(+, view(GQweights, idx), view(GQnodes, idx)) do weight, node -# weight * f(node) -# end -# end -# end -# mapreduce(fetch, +, tasks) -# end -# end + betterSum(GQweights .* BroadcastArray(genFγIntegrand(γ, u), GQnodes)) +end + function F0(u::T) where {T} if u < getAtolVal(T) @@ -79,14 +58,13 @@ function Fγ(γ::Int, u::T) where {T} end end -function F₀toFγ(γ::Int, u::T) where {T} - res = Array{T}(undef, γ+1) - res[begin] = F0(u) - γ > 0 && (res[end] = Fγ(γ, u)) +function F₀toFγ(γ::Int, u::T, resHolder::Vector{T}=Array{T}(undef, γ+1)) where {T} + resHolder[begin] = F0(u) + γ > 0 && (resHolder[γ+1] = Fγ(γ, u)) for i in γ:-1:2 - @inbounds res[i] = (expm1(-u) + 2u*res[i+1] + 1) / (2i - 1) + @inbounds resHolder[i] = (expm1(-u) + 2u*resHolder[i+1] + 1) / (2i - 1) end - res + resHolder end @@ -179,6 +157,8 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T A = T(0.0) i₁, j₁, k₁ = ijk₁ i₂, j₂, k₂ = ijk₂ + # ijkSum = sum(ijk₁) + sum(ijk₂) + # Fγss = [F₀toFγ(γ, β) for γ in 0:ijkSum] for l₁ in 0:(i₁÷2), m₁ in 0:(j₁÷2), n₁ in 0:(k₁÷2), l₂ in 0:(i₂÷2), m₂ in 0:(j₂÷2), n₂ in 0:(k₂÷2) @@ -194,6 +174,7 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq₁+opq₂) μsum = sum(μv) Fγs = F₀toFγ(μsum, β) + # Fγs = Fγss[μsum+1] core1s = genIntTerm1.(ΔR₁R₂, lmn₁, opq₁, lmn₂, opq₂, ijk₁, α₁, ijk₂, α₂) for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) diff --git a/src/Tools.jl b/src/Tools.jl index e0f8fec9..4b19ba19 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -1007,4 +1007,24 @@ function keepOnly!(a::AbstractArray, idx::Int) e = a[idx] deleteat!(a, (firstindex(a)+1):lastindex(a)) a[] = e +end + + +function betterSum(x) # Improved Kahan–Babuška algorithm fro array summation + xFirxt = first(x) + s = xFirxt .- xFirxt + c = xFirxt .- xFirxt + + for xi in x + t = s .+ xi + c += map(s, t, xi) do si, ti, xij + if abs(si) >= abs(xij) + (si - ti) + xij + else + (xij - ti) + si + end + end + s = t + end + s .+ c end \ No newline at end of file From 71e3275e92768ee7026ea311f1a0e91c40b0224f Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Tue, 2 Jan 2024 18:12:24 -0500 Subject: [PATCH 08/22] Performance optimization. --- src/Integrals/Core.jl | 48 ++++++++++++++++++++++++------------------- 1 file changed, 27 insertions(+), 21 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 9385fdcb..ff9c2480 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -72,17 +72,17 @@ function genIntOverlapCore(Δx::T, i₁::Int, α₁::T, i₂::Int, α₂::T) where {T} res = T(0.0) + getRange = ifelse( iszero(Δx), + Ω::Int -> ifelse(iseven(Ω), begin halfΩ=Ω÷2; halfΩ:halfΩ end, 1:0), + Ω::Int -> 0:(Ω÷2) ) for l₁ in 0:(i₁÷2), l₂ in 0:(i₂÷2) Ω = muladd(-2, l₁ + l₂, i₁ + i₂) - halfΩ = Ω÷2 - oRange = 0:halfΩ - Δx == 0.0 && (iseven(Ω) ? (oRange = halfΩ:halfΩ) : continue) - for o in oRange + for o in getRange(Ω) res += Δx^(Ω-2o) * ( α₁^(i₂ - l₁ - 2l₂ - o) / (factorial(l₂) * factorial(i₁-2l₁)) ) * ( α₂^(i₁ - l₂ - 2l₁ - o) / (factorial(l₁) * factorial(i₂-2l₂)) ) * - T( (-1)^o * factorial(Ω) / - (4^(l₁+ l₂ + o) * factorial(o) * factorial(Ω-2o)) ) * + ( T(-1)^o * factorial(Ω) / + (( 1 << 2(l₁+ l₂ + o) ) * factorial(o ) * factorial(Ω -2o )) ) * (α₁ + α₂)^muladd(2, (l₁ + l₂), o) end end @@ -92,18 +92,18 @@ end function ∫overlapCore(::Val{3}, ΔR::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, - ijk₂::NTuple{3, Int}, α₂::T) where {T} + ijk₂::NTuple{3, Int}, α₂::T, coeff::T=T(1)) where {T} any(n -> n<0, (ijk₁..., ijk₂...)) && (return T(0.0)) α = α₁ + α₂ res = T(1) for (i₁, i₂, ΔRᵢ) in zip(ijk₁, ijk₂, ΔR) int = genIntOverlapCore(ΔRᵢ, i₁, α₁, i₂, α₂) - iszero(int) && (return T(0)) + iszero(int) && (return T(0.0)) res *= (-1)^(i₁) * factorial(i₁) * factorial(i₂) * α^(-i₁-i₂) * int end - res *= sqrt((π/α)^3) * exp(-α₁ / α * α₂* sum(abs2, ΔR)) - res + res *= sqrt((π/α)^3) * exp(-α₁ / α * α₂ * sum(abs2, ΔR)) + res * coeff end ∫overlapCore(::Val{3}, @@ -120,10 +120,15 @@ function ∫elecKineticCore(::Val{3}, ΔR = R₁ .- R₂ shifts = ((1,0,0), (0,1,0), (0,0,1)) mapreduce(+, ijk₁, ijk₂, shifts) do 𝑙₁c, 𝑙₂c, Δ𝑙 - ∫overlapCore(Val(3), ΔR, map(-, ijk₁, Δ𝑙), α₁, map(-, ijk₂, Δ𝑙), α₂) * 𝑙₁c*𝑙₂c/2 + - ∫overlapCore(Val(3), ΔR, map(+, ijk₁, Δ𝑙), α₁, map(+, ijk₂, Δ𝑙), α₂) * 2α₁*α₂ - - ∫overlapCore(Val(3), ΔR, map(+, ijk₁, Δ𝑙), α₁, map(-, ijk₂, Δ𝑙), α₂) * α₁*𝑙₂c - - ∫overlapCore(Val(3), ΔR, map(-, ijk₁, Δ𝑙), α₁, map(+, ijk₂, Δ𝑙), α₂) * α₂*𝑙₁c + Δijk1 = map(-, ijk₁, Δ𝑙) + Δijk2 = map(-, ijk₂, Δ𝑙) + Δijk3 = map(+, ijk₁, Δ𝑙) + Δijk4 = map(+, ijk₂, Δ𝑙) + int1 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk2, α₂, T(𝑙₁c*𝑙₂c/2)) + int2 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk4, α₂, 2α₁*α₂) + int3 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk2, α₂, α₁*𝑙₂c) + int4 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk4, α₂, α₂*𝑙₁c) + int1 + int2 - int3 - int4 end end @@ -134,15 +139,17 @@ function genIntTerm1(Δx::T1, i₁::T2, α₁::T1, i₂::T2, α₂::T1) where {T1, T2<:Integer} (r::T2) -> - ( Δx^muladd(-2, r, o₁+o₂) / (factorial(r ) * (factorial∘muladd)(-2, r, o₁+o₂)) ) * + ( Δx^muladd(-2, r, o₁+o₂) / (factorial(r) * (factorial∘muladd)(-2, r, o₁+o₂)) ) * ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * ( α₂^(o₁-l₂- r) / (factorial(l₂) * factorial(i₂-2l₂-o₂)) ) * - T1( (-1)^(o₂+r) * factorial(o₁+o₂) / (4^(l₁+l₂+r) * factorial(o₁) * factorial(o₂)) ) + T1( (-1)^(o₂+r) * factorial(o₁+o₂) / + ((1 << 2(l₁+l₂+r)) * factorial(o₁) * factorial(o₂)) ) end function genIntTerm2core(Δx::T1, μ::T2) where {T1, T2<:Integer} (u::T2) -> - Δx^(μ-2u) * T1( (-1)^u * factorial(μ) / (4^u * factorial(u) * factorial(μ-2u)) ) + Δx^(μ-2u) * T1( (-1)^u * factorial(μ) / + ((1 << 2u) * factorial(u) * factorial(μ-2u)) ) end function genIntTerm2(Δx::T1, α::T1, o₁::T2, o₂::T2, μ::T2, r::T2) where {T1, T2<:Integer} @@ -157,8 +164,8 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T A = T(0.0) i₁, j₁, k₁ = ijk₁ i₂, j₂, k₂ = ijk₂ - # ijkSum = sum(ijk₁) + sum(ijk₂) - # Fγss = [F₀toFγ(γ, β) for γ in 0:ijkSum] + ijkSum = sum(ijk₁) + sum(ijk₂) + Fγss = [F₀toFγ(γ, β) for γ in 0:ijkSum] for l₁ in 0:(i₁÷2), m₁ in 0:(j₁÷2), n₁ in 0:(k₁÷2), l₂ in 0:(i₂÷2), m₂ in 0:(j₂÷2), n₂ in 0:(k₂÷2) @@ -173,8 +180,7 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq₁+opq₂) μsum = sum(μv) - Fγs = F₀toFγ(μsum, β) - # Fγs = Fγss[μsum+1] + Fγs = @inbounds Fγss[μsum+1] core1s = genIntTerm1.(ΔR₁R₂, lmn₁, opq₁, lmn₂, opq₂, ijk₁, α₁, ijk₂, α₂) for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) From bd5457c79c6edb542c23a0a0a5de2729fc6df4bb Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Tue, 2 Jan 2024 23:57:53 -0500 Subject: [PATCH 09/22] Divided `genIntOverlapCore` into sub functions for perf opt. --- src/Integrals/Core.jl | 176 +++++++++++++++++++++++++++++------------- 1 file changed, 121 insertions(+), 55 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index ff9c2480..815ad611 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -37,7 +37,7 @@ end function F0(u::T) where {T} if u < getAtolVal(T) - T(1) + T(1.0) else ur = sqrt(u) T(πvals[0.5]) * erf(ur) / (2ur) @@ -68,37 +68,95 @@ function F₀toFγ(γ::Int, u::T, resHolder::Vector{T}=Array{T}(undef, γ+1)) wh end -function genIntOverlapCore(Δx::T, - i₁::Int, α₁::T, - i₂::Int, α₂::T) where {T} +function genIntOverlapTermAA(getRangeFunc::F, Δx::T, i::Int, α::T) where {F, T} + f = let getRange=getRangeFunc, Δx=Δx, i=i, α=α + function (l₁::Int, l₂::Int) + res = T(0.0) + Ω = muladd(-2, l₁ + l₂, 2i) + for o in getRange(Ω) + res += Δx^(Ω-2o) * + ( α^(2i - 3l₁ - 3l₂ - 2o) / (factorial(l₂) * factorial(i-2l₁) * + factorial(l₁) * factorial(i-2l₂)) ) * + ( T(-1)^o * factorial(Ω) / + ((1 << 2(l₁ + l₂ + o)) * factorial(o ) * factorial(Ω-2o )) ) * + (2α)^muladd(2, (l₁ + l₂), o) + end + res + end + end + f +end + +function genIntOverlapTermAB(getRangeFunc::F, Δx::T, i₁::Int, α₁::T, + i₂::Int, α₂::T) where {F, T} + f = let getRange=getRangeFunc, Δx=Δx, i₁=i₁, α₁=α₁, i₂=i₂, α₂=α₂ + function (l₁::Int, l₂::Int) + res = T(0.0) + Ω = muladd(-2, l₁ + l₂, i₁ + i₂) + for o in getRange(Ω) + res += Δx^(Ω-2o) * + ( α₁^(i₂ - l₁ - 2l₂ - o) / (factorial(l₂) * factorial(i₁-2l₁)) ) * + ( α₂^(i₁ - l₂ - 2l₁ - o) / (factorial(l₁) * factorial(i₂-2l₂)) ) * + ( T(-1)^o * factorial(Ω) / + ((1 << 2(l₁ + l₂ + o)) * factorial(o ) * factorial(Ω -2o )) ) * + (α₁ + α₂)^muladd(2, (l₁ + l₂), o) + end + res + end + end + f +end + +function genIntOverlapCore11(Δx::T, i::Int, α::T, + sameCen::Bool=false) where {T} + getRange = ifelse( sameCen || iszero(Δx), + Ω::Int -> ifelse(iseven(Ω), begin ΩHalf=Ω÷2; ΩHalf:ΩHalf end, 1:0), + Ω::Int -> 0:(Ω÷2) ) + f = genIntOverlapTermAA(getRange, Δx, i, α) res = T(0.0) - getRange = ifelse( iszero(Δx), - Ω::Int -> ifelse(iseven(Ω), begin halfΩ=Ω÷2; halfΩ:halfΩ end, 1:0), + iHalf = i÷2 + for l₁ in 0:iHalf, l₂ in l₁:iHalf + res += f(l₁, l₂) * (1 + (l₁!=l₂)) + end + res +end + + +function genIntOverlapCore12(Δx::T, i₁::Int, α₁::T, + i₂::Int, α₂::T, + sameCen::Bool=false) where {T} + getRange = ifelse( sameCen || iszero(Δx), + Ω::Int -> ifelse(iseven(Ω), begin ΩHalf=Ω÷2; ΩHalf:ΩHalf end, 1:0), Ω::Int -> 0:(Ω÷2) ) + f = genIntOverlapTermAB(getRange, Δx, i₁, α₁, i₂, α₂) + res = T(0.0) for l₁ in 0:(i₁÷2), l₂ in 0:(i₂÷2) - Ω = muladd(-2, l₁ + l₂, i₁ + i₂) - for o in getRange(Ω) - res += Δx^(Ω-2o) * - ( α₁^(i₂ - l₁ - 2l₂ - o) / (factorial(l₂) * factorial(i₁-2l₁)) ) * - ( α₂^(i₁ - l₂ - 2l₁ - o) / (factorial(l₁) * factorial(i₂-2l₂)) ) * - ( T(-1)^o * factorial(Ω) / - (( 1 << 2(l₁+ l₂ + o) ) * factorial(o ) * factorial(Ω -2o )) ) * - (α₁ + α₂)^muladd(2, (l₁ + l₂), o) - end + res += f(l₁, l₂) end res end +function genIntOverlapCore(Δx::T, i₁::Int, α₁::T, + i₂::Int, α₂::T, + (sameCenAngXpn)::NTuple{3, Bool}=(false, false, false)) where {T} + if (sameCenAngXpn[2] || i₁==i₂) && (sameCenAngXpn[3] || α₁==α₂) + genIntOverlapCore11(Δx, i₁, α₁, sameCenAngXpn[1]) + else + genIntOverlapCore12(Δx, i₁, α₁, i₂, α₂, sameCenAngXpn[1]) + end +end + function ∫overlapCore(::Val{3}, ΔR::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, - ijk₂::NTuple{3, Int}, α₂::T, coeff::T=T(1)) where {T} + ijk₂::NTuple{3, Int}, α₂::T, + sameCenAngXpn::NTuple{3, Bool}=(false, false, false), + coeff::T=T(1.0)) where {T} any(n -> n<0, (ijk₁..., ijk₂...)) && (return T(0.0)) - α = α₁ + α₂ - res = T(1) + res = T(1.0) for (i₁, i₂, ΔRᵢ) in zip(ijk₁, ijk₂, ΔR) - int = genIntOverlapCore(ΔRᵢ, i₁, α₁, i₂, α₂) + int = genIntOverlapCore(ΔRᵢ, i₁, α₁, i₂, α₂, sameCenAngXpn) iszero(int) && (return T(0.0)) res *= (-1)^(i₁) * factorial(i₁) * factorial(i₂) * α^(-i₁-i₂) * int end @@ -109,35 +167,51 @@ end ∫overlapCore(::Val{3}, R₁::NTuple{3, T}, R₂::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, - ijk₂::NTuple{3, Int}, α₂::T) where {T} = -∫overlapCore(Val(3), R₁.-R₂, ijk₁, α₁, ijk₂, α₂) + ijk₂::NTuple{3, Int}, α₂::T, + sameCenAngXpn::NTuple{3, Bool}=(false, false, false), + coeff::T=T(1.0)) where {T} = +∫overlapCore(Val(3), R₁.-R₂, ijk₁, α₁, ijk₂, α₂, sameCenAngXpn, coeff) function ∫elecKineticCore(::Val{3}, R₁::NTuple{3, T}, R₂::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, - ijk₂::NTuple{3, Int}, α₂::T) where {T} + ijk₂::NTuple{3, Int}, α₂::T, + (blC, _, blX)::NTuple{3, Bool}=(false, false, false)) where {T} + sameCenAngXpn = (blC, false, blX) ΔR = R₁ .- R₂ shifts = ((1,0,0), (0,1,0), (0,0,1)) - mapreduce(+, ijk₁, ijk₂, shifts) do 𝑙₁c, 𝑙₂c, Δ𝑙 + mapreduce(+, ijk₁, ijk₂, shifts) do i₁, i₂, Δ𝑙 Δijk1 = map(-, ijk₁, Δ𝑙) Δijk2 = map(-, ijk₂, Δ𝑙) Δijk3 = map(+, ijk₁, Δ𝑙) Δijk4 = map(+, ijk₂, Δ𝑙) - int1 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk2, α₂, T(𝑙₁c*𝑙₂c/2)) - int2 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk4, α₂, 2α₁*α₂) - int3 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk2, α₂, α₁*𝑙₂c) - int4 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk4, α₂, α₂*𝑙₁c) + int1 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk2, α₂, sameCenAngXpn, T(i₁*i₂/2)) + int2 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk4, α₂, sameCenAngXpn, 2α₁*α₂ ) + int3 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk2, α₂, sameCenAngXpn, α₁*i₂ ) + int4 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk4, α₂, sameCenAngXpn, α₂*i₁ ) int1 + int2 - int3 - int4 end end -function genIntTerm1(Δx::T1, - l₁::T2, o₁::T2, - l₂::T2, o₂::T2, - i₁::T2, α₁::T1, - i₂::T2, α₂::T1) where {T1, T2<:Integer} +# function genIntTerm1AA(Δx::T1, +# l₁::T2, o₁::T2, +# l₂::T2, o₂::T2, +# i::T2, α::T1) where {T1, T2<:Integer} +# (r::T2) -> +# ( Δx^muladd(-2, r, o₁+o₂) / (factorial(r) * (factorial∘muladd)(-2, r, o₁+o₂)) ) * +# ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * +# ( α₂^(o₁-l₂- r) / (factorial(l₂) * factorial(i₂-2l₂-o₂)) ) * +# T1( (-1)^(o₂+r) * factorial(o₁+o₂) / +# ((1 << 2(l₁+l₂+r)) * factorial(o₁) * factorial(o₂)) ) +# end + +function genIntTerm1AB(Δx::T1, + l₁::T2, o₁::T2, + l₂::T2, o₂::T2, + i₁::T2, α₁::T1, + i₂::T2, α₂::T1) where {T1, T2<:Integer} (r::T2) -> ( Δx^muladd(-2, r, o₁+o₂) / (factorial(r) * (factorial∘muladd)(-2, r, o₁+o₂)) ) * ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * @@ -160,7 +234,8 @@ end function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T, ijk₁::NTuple{3, Int}, α₁::T, - ijk₂::NTuple{3, Int}, α₂::T) where {T} + ijk₂::NTuple{3, Int}, α₂::T, + sameCenAngXpn::NTuple{3, Bool}=(false, false, false)) where {T} A = T(0.0) i₁, j₁, k₁ = ijk₁ i₂, j₂, k₂ = ijk₂ @@ -181,7 +256,7 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq₁+opq₂) μsum = sum(μv) Fγs = @inbounds Fγss[μsum+1] - core1s = genIntTerm1.(ΔR₁R₂, lmn₁, opq₁, lmn₂, opq₂, ijk₁, α₁, ijk₂, α₂) + core1s = genIntTerm1AB.(ΔR₁R₂, lmn₁, opq₁, lmn₂, opq₂, ijk₁, α₁, ijk₂, α₂) for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) @@ -206,20 +281,14 @@ function ∫nucAttractionCore(::Val{3}, Z₀::Int, R₀::NTuple{3, T}, R₁::NTuple{3, T}, R₂::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, - ijk₂::NTuple{3, Int}, α₂::T) where {T} - if α₁ == α₂ - α = 2α₁ - R = @. (R₁ + R₂) / 2 - flag = true - else - α = α₁ + α₂ - R = @. (α₁*R₁ + α₂*R₂) / α - flag = false - end + ijk₂::NTuple{3, Int}, α₂::T, + sameCenAngXpn::NTuple{3, Bool}=(false, false, false)) where {T} + α = α₁ + α₂ + R = @. (α₁*R₁ + α₂*R₂) / α ΔRR₀ = R .- R₀ ΔR₁R₂ = R₁ .- R₂ β = α * sum(abs2, ΔRR₀) - genIntNucAttCore(ΔRR₀, ΔR₁R₂, β, ijk₁, α₁, ijk₂, α₂) * + genIntNucAttCore(ΔRR₀, ΔR₁R₂, β, ijk₁, α₁, ijk₂, α₂, sameCenAngXpn) * (π / α) * exp(-α₁ / α * α₂ * sum(abs2, ΔR₁R₂)) * ( -Z₀ * (-1)^sum(ijk₁ .+ ijk₂) * mapMapReduce(ijk₁, factorial) * mapMapReduce(ijk₂, factorial) ) @@ -231,7 +300,7 @@ function genIntTerm3(Δx::T1, i₁::T2, α₁::T1, i₂::T2, α₂::T1) where {T1, T2<:Integer} (r::T2) -> - genIntTerm1(Δx, l₁, o₁, l₂, o₂, i₁, α₁, i₂, α₂)(r) * (α₁+α₂)^muladd(2, l₁+l₂, r) + genIntTerm1AB(Δx, l₁, o₁, l₂, o₂, i₁, α₁, i₂, α₂)(r) * (α₁+α₂)^muladd(2, l₁+l₂, r) end function genIntTerm4(Δx::T1, η::T1, μ::T2) where {T1, T2<:Integer} @@ -458,13 +527,8 @@ end diFoldCount(i::T, j::T) where {T} = ifelse(i==j, 1, 2) -function octaFoldCount(i::T, j::T, k::T, l::T) where {T} - m = 0 - i != j && (m += 1) - k != l && (m += 1) - (i != k || j != l) && (m += 1) - 1 << m -end +octaFoldCount(i::T, j::T, k::T, l::T) where {T} = +1 << ( (i != j) + (k != l) + (i != k || j != l) ) function getIntCore11!(n::Int, @@ -566,9 +630,11 @@ function getOneBodyInt(::Type{T}, ::Val{D}, ∫1e::F, @nospecialize(optPosArgs:: {T, D, F<:Function} (R₁, ijk₁, ps₁, 𝑙₁), (R₂, ijk₂, ps₂, 𝑙₂) = reformatIntData1(iBl, bfs) 𝑙₁==𝑙₂==0 || isIntZero(F, optPosArgs, R₁,R₂, ijk₁,ijk₂) && (return T(0.0)) - uniquePairs, uPairCoeffs = getOneBodyUniquePairs(R₁==R₂ && ijk₁==ijk₂, ps₁, ps₂) + sameCenAngXpn = (R₁==R₂, ijk₁==ijk₂, false) # There might be a way to get the 3rd entry. + uniquePairs, uPairCoeffs = getOneBodyUniquePairs(sameCenAngXpn[1] && sameCenAngXpn[2], + ps₁, ps₂) mapreduce(+, uniquePairs, uPairCoeffs) do x, y - ∫1e(Val(D), optPosArgs..., R₁, R₂, ijk₁, x[1], ijk₂, x[2])::T * y + ∫1e(Val(D), optPosArgs..., R₁, R₂, ijk₁, x[1], ijk₂, x[2], sameCenAngXpn)::T * y end end From 7ae88d6961e47e84265cc6302af969d82d4cb510 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Wed, 3 Jan 2024 04:56:01 -0500 Subject: [PATCH 10/22] Removed ineefective opt and added `sameCenAng`. --- src/Integrals/Core.jl | 165 ++++++++++++++++++++---------------------- 1 file changed, 79 insertions(+), 86 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 815ad611..66a869b9 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -68,27 +68,8 @@ function F₀toFγ(γ::Int, u::T, resHolder::Vector{T}=Array{T}(undef, γ+1)) wh end -function genIntOverlapTermAA(getRangeFunc::F, Δx::T, i::Int, α::T) where {F, T} - f = let getRange=getRangeFunc, Δx=Δx, i=i, α=α - function (l₁::Int, l₂::Int) - res = T(0.0) - Ω = muladd(-2, l₁ + l₂, 2i) - for o in getRange(Ω) - res += Δx^(Ω-2o) * - ( α^(2i - 3l₁ - 3l₂ - 2o) / (factorial(l₂) * factorial(i-2l₁) * - factorial(l₁) * factorial(i-2l₂)) ) * - ( T(-1)^o * factorial(Ω) / - ((1 << 2(l₁ + l₂ + o)) * factorial(o ) * factorial(Ω-2o )) ) * - (2α)^muladd(2, (l₁ + l₂), o) - end - res - end - end - f -end - -function genIntOverlapTermAB(getRangeFunc::F, Δx::T, i₁::Int, α₁::T, - i₂::Int, α₂::T) where {F, T} +function genIntOverlapTerm(getRangeFunc::F, Δx::T, i₁::Int, α₁::T, + i₂::Int, α₂::T) where {F, T} f = let getRange=getRangeFunc, Δx=Δx, i₁=i₁, α₁=α₁, i₂=i₂, α₂=α₂ function (l₁::Int, l₂::Int) res = T(0.0) @@ -107,28 +88,13 @@ function genIntOverlapTermAB(getRangeFunc::F, Δx::T, i₁::Int, α₁::T, f end -function genIntOverlapCore11(Δx::T, i::Int, α::T, - sameCen::Bool=false) where {T} - getRange = ifelse( sameCen || iszero(Δx), - Ω::Int -> ifelse(iseven(Ω), begin ΩHalf=Ω÷2; ΩHalf:ΩHalf end, 1:0), - Ω::Int -> 0:(Ω÷2) ) - f = genIntOverlapTermAA(getRange, Δx, i, α) - res = T(0.0) - iHalf = i÷2 - for l₁ in 0:iHalf, l₂ in l₁:iHalf - res += f(l₁, l₂) * (1 + (l₁!=l₂)) - end - res -end - - -function genIntOverlapCore12(Δx::T, i₁::Int, α₁::T, - i₂::Int, α₂::T, - sameCen::Bool=false) where {T} +function genIntOverlapCore(Δx::T, i₁::Int, α₁::T, + i₂::Int, α₂::T, + sameCen::Bool=false) where {T} getRange = ifelse( sameCen || iszero(Δx), Ω::Int -> ifelse(iseven(Ω), begin ΩHalf=Ω÷2; ΩHalf:ΩHalf end, 1:0), Ω::Int -> 0:(Ω÷2) ) - f = genIntOverlapTermAB(getRange, Δx, i₁, α₁, i₂, α₂) + f = genIntOverlapTerm(getRange, Δx, i₁, α₁, i₂, α₂) res = T(0.0) for l₁ in 0:(i₁÷2), l₂ in 0:(i₂÷2) res += f(l₁, l₂) @@ -136,27 +102,17 @@ function genIntOverlapCore12(Δx::T, i₁::Int, α₁::T, res end -function genIntOverlapCore(Δx::T, i₁::Int, α₁::T, - i₂::Int, α₂::T, - (sameCenAngXpn)::NTuple{3, Bool}=(false, false, false)) where {T} - if (sameCenAngXpn[2] || i₁==i₂) && (sameCenAngXpn[3] || α₁==α₂) - genIntOverlapCore11(Δx, i₁, α₁, sameCenAngXpn[1]) - else - genIntOverlapCore12(Δx, i₁, α₁, i₂, α₂, sameCenAngXpn[1]) - end -end - function ∫overlapCore(::Val{3}, ΔR::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - sameCenAngXpn::NTuple{3, Bool}=(false, false, false), + sameCenAng::NTuple{2, Bool}=(false, false), coeff::T=T(1.0)) where {T} any(n -> n<0, (ijk₁..., ijk₂...)) && (return T(0.0)) α = α₁ + α₂ res = T(1.0) for (i₁, i₂, ΔRᵢ) in zip(ijk₁, ijk₂, ΔR) - int = genIntOverlapCore(ΔRᵢ, i₁, α₁, i₂, α₂, sameCenAngXpn) + int = genIntOverlapCore(ΔRᵢ, i₁, α₁, i₂, α₂, sameCenAng[begin]) iszero(int) && (return T(0.0)) res *= (-1)^(i₁) * factorial(i₁) * factorial(i₂) * α^(-i₁-i₂) * int end @@ -168,17 +124,17 @@ end R₁::NTuple{3, T}, R₂::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - sameCenAngXpn::NTuple{3, Bool}=(false, false, false), + sameCenAng::NTuple{2, Bool}=(false, false), coeff::T=T(1.0)) where {T} = -∫overlapCore(Val(3), R₁.-R₂, ijk₁, α₁, ijk₂, α₂, sameCenAngXpn, coeff) +∫overlapCore(Val(3), R₁.-R₂, ijk₁, α₁, ijk₂, α₂, sameCenAng, coeff) function ∫elecKineticCore(::Val{3}, R₁::NTuple{3, T}, R₂::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - (blC, _, blX)::NTuple{3, Bool}=(false, false, false)) where {T} - sameCenAngXpn = (blC, false, blX) + sameCenAng::NTuple{2, Bool}=(false, false)) where {T} + sameCenAng = (sameCenAng[begin], false) ΔR = R₁ .- R₂ shifts = ((1,0,0), (0,1,0), (0,0,1)) mapreduce(+, ijk₁, ijk₂, shifts) do i₁, i₂, Δ𝑙 @@ -186,32 +142,20 @@ function ∫elecKineticCore(::Val{3}, Δijk2 = map(-, ijk₂, Δ𝑙) Δijk3 = map(+, ijk₁, Δ𝑙) Δijk4 = map(+, ijk₂, Δ𝑙) - int1 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk2, α₂, sameCenAngXpn, T(i₁*i₂/2)) - int2 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk4, α₂, sameCenAngXpn, 2α₁*α₂ ) - int3 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk2, α₂, sameCenAngXpn, α₁*i₂ ) - int4 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk4, α₂, sameCenAngXpn, α₂*i₁ ) + int1 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk2, α₂, sameCenAng, T(i₁*i₂/2)) + int2 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk4, α₂, sameCenAng, 2α₁*α₂ ) + int3 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk2, α₂, sameCenAng, α₁*i₂ ) + int4 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk4, α₂, sameCenAng, α₂*i₁ ) int1 + int2 - int3 - int4 end end -# function genIntTerm1AA(Δx::T1, -# l₁::T2, o₁::T2, -# l₂::T2, o₂::T2, -# i::T2, α::T1) where {T1, T2<:Integer} -# (r::T2) -> -# ( Δx^muladd(-2, r, o₁+o₂) / (factorial(r) * (factorial∘muladd)(-2, r, o₁+o₂)) ) * -# ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * -# ( α₂^(o₁-l₂- r) / (factorial(l₂) * factorial(i₂-2l₂-o₂)) ) * -# T1( (-1)^(o₂+r) * factorial(o₁+o₂) / -# ((1 << 2(l₁+l₂+r)) * factorial(o₁) * factorial(o₂)) ) -# end - -function genIntTerm1AB(Δx::T1, - l₁::T2, o₁::T2, - l₂::T2, o₂::T2, - i₁::T2, α₁::T1, - i₂::T2, α₂::T1) where {T1, T2<:Integer} +function genIntTerm1(Δx::T1, + l₁::T2, o₁::T2, + l₂::T2, o₂::T2, + i₁::T2, α₁::T1, + i₂::T2, α₂::T1) where {T1, T2<:Integer} (r::T2) -> ( Δx^muladd(-2, r, o₁+o₂) / (factorial(r) * (factorial∘muladd)(-2, r, o₁+o₂)) ) * ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * @@ -235,7 +179,7 @@ end function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T, ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - sameCenAngXpn::NTuple{3, Bool}=(false, false, false)) where {T} + sameCenAng::NTuple{2, Bool}=(false, false)) where {T} A = T(0.0) i₁, j₁, k₁ = ijk₁ i₂, j₂, k₂ = ijk₂ @@ -256,7 +200,7 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq₁+opq₂) μsum = sum(μv) Fγs = @inbounds Fγss[μsum+1] - core1s = genIntTerm1AB.(ΔR₁R₂, lmn₁, opq₁, lmn₂, opq₂, ijk₁, α₁, ijk₂, α₂) + core1s = genIntTerm1.(ΔR₁R₂, lmn₁, opq₁, lmn₂, opq₂, ijk₁, α₁, ijk₂, α₂) for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) @@ -276,19 +220,69 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T A end +# function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T, +# ijk₁::NTuple{3, Int}, α₁::T, +# ijk₂::NTuple{3, Int}, α₂::T, +# sameCenAng::NTuple{2, Bool}=(false, false)) where {T} +# A = T(0.0) +# i₁, j₁, k₁ = ijk₁ +# i₂, j₂, k₂ = ijk₂ +# ijkSum = sum(ijk₁) + sum(ijk₂) +# Fγss = [F₀toFγ(γ, β) for γ in 0:ijkSum] +# for l₁ in 0:(i₁÷2), m₁ in 0:(j₁÷2), n₁ in 0:(k₁÷2), +# l₂ in 0:(i₂÷2), m₂ in 0:(j₂÷2), n₂ in 0:(k₂÷2) + +# lmn₁ = (l₁, m₁, n₁) +# lmn₂ = (l₂, m₂, n₂) + +# res = T(0) + +# for o₁ in 0:(i₁-2l₁), p₁ in 0:(j₁-2m₁), q₁ in 0:(k₁-2n₁), +# o₂ in 0:(i₂-2l₂), p₂ in 0:(j₂-2m₂), q₂ in 0:(k₂-2n₂) + +# opq₁ = (o₁, p₁, q₁) +# opq₂ = (o₂, p₂, q₂) + +# μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq₁+opq₂) +# μsum = sum(μv) +# Fγs = @inbounds Fγss[μsum+1] +# core1s = genIntTerm1.(ΔR₁R₂, lmn₁, opq₁, lmn₂, opq₂, ijk₁, α₁, ijk₂, α₂) + +# for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) + +# rst = (r, s, t) +# tmp = T(0.0) +# core2s = genIntTerm2.(ΔRR₀, α₁+α₂, opq₁, opq₂, μv, rst) + +# for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) +# γ = μsum - u - v - w +# @inbounds tmp += mapMapReduce((u, v, w), core2s) * 2Fγs[γ+1] +# end +# rr1 = mapMapReduce(rst, core1s) +# rr2 = tmp +# println(((rr1, (o₁+o₂)÷2, (p₁+p₂)÷2, (q₁+q₂)÷2), (rr2, μˣ, μʸ, μᶻ)), ",") +# res += rr1 * rr2 +# end +# end +# A += res +# end +# println() +# A +# end + function ∫nucAttractionCore(::Val{3}, Z₀::Int, R₀::NTuple{3, T}, R₁::NTuple{3, T}, R₂::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - sameCenAngXpn::NTuple{3, Bool}=(false, false, false)) where {T} + sameCenAng::NTuple{2, Bool}=(false, false)) where {T} α = α₁ + α₂ R = @. (α₁*R₁ + α₂*R₂) / α ΔRR₀ = R .- R₀ ΔR₁R₂ = R₁ .- R₂ β = α * sum(abs2, ΔRR₀) - genIntNucAttCore(ΔRR₀, ΔR₁R₂, β, ijk₁, α₁, ijk₂, α₂, sameCenAngXpn) * + genIntNucAttCore(ΔRR₀, ΔR₁R₂, β, ijk₁, α₁, ijk₂, α₂, sameCenAng) * (π / α) * exp(-α₁ / α * α₂ * sum(abs2, ΔR₁R₂)) * ( -Z₀ * (-1)^sum(ijk₁ .+ ijk₂) * mapMapReduce(ijk₁, factorial) * mapMapReduce(ijk₂, factorial) ) @@ -300,7 +294,7 @@ function genIntTerm3(Δx::T1, i₁::T2, α₁::T1, i₂::T2, α₂::T1) where {T1, T2<:Integer} (r::T2) -> - genIntTerm1AB(Δx, l₁, o₁, l₂, o₂, i₁, α₁, i₂, α₂)(r) * (α₁+α₂)^muladd(2, l₁+l₂, r) + genIntTerm1(Δx, l₁, o₁, l₂, o₂, i₁, α₁, i₂, α₂)(r) * (α₁+α₂)^muladd(2, l₁+l₂, r) end function genIntTerm4(Δx::T1, η::T1, μ::T2) where {T1, T2<:Integer} @@ -630,11 +624,10 @@ function getOneBodyInt(::Type{T}, ::Val{D}, ∫1e::F, @nospecialize(optPosArgs:: {T, D, F<:Function} (R₁, ijk₁, ps₁, 𝑙₁), (R₂, ijk₂, ps₂, 𝑙₂) = reformatIntData1(iBl, bfs) 𝑙₁==𝑙₂==0 || isIntZero(F, optPosArgs, R₁,R₂, ijk₁,ijk₂) && (return T(0.0)) - sameCenAngXpn = (R₁==R₂, ijk₁==ijk₂, false) # There might be a way to get the 3rd entry. - uniquePairs, uPairCoeffs = getOneBodyUniquePairs(sameCenAngXpn[1] && sameCenAngXpn[2], - ps₁, ps₂) + sameCenAng = (R₁==R₂, ijk₁==ijk₂) + uniquePairs, uPairCoeffs = getOneBodyUniquePairs(prod(sameCenAng), ps₁, ps₂) mapreduce(+, uniquePairs, uPairCoeffs) do x, y - ∫1e(Val(D), optPosArgs..., R₁, R₂, ijk₁, x[1], ijk₂, x[2], sameCenAngXpn)::T * y + ∫1e(Val(D), optPosArgs..., R₁, R₂, ijk₁, x[1], ijk₂, x[2], sameCenAng)::T * y end end From 35c627017586da9b43a5e57734f4c334203b6234 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Sun, 7 Jan 2024 23:33:47 -0500 Subject: [PATCH 11/22] Performance optimization. --- src/CI.jl | 2 +- src/Integrals/Core.jl | 182 +++++++++++++++++++++++++++--------------- 2 files changed, 120 insertions(+), 64 deletions(-) diff --git a/src/CI.jl b/src/CI.jl index 5bfb4c90..52361d9e 100644 --- a/src/CI.jl +++ b/src/CI.jl @@ -19,7 +19,7 @@ struct CSF{𝑠, 𝑚ₛ, ON} <: RCI{𝑚ₛ, ON} function RSD(spatialOccu::NTuple{ON, Bool}, basis::NTuple{ON, <:GTBasisFuncs{T, D, 1}}) - for + # for new{Nup-Ndn, ON}(Nup+Ndn, spinUpOccu, spinDnOccu, objectid(basis)) end end diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 66a869b9..ff4b2e1d 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -12,9 +12,9 @@ using LazyArrays: BroadcastArray ## [DOI] 10.48550/arXiv.2007.12057 function genFγIntegrand(γ::Int, u::T) where {T} - f = let γLoc=γ, uLoc=u + f = let γ=γ, u=u @inline function (x::T) # @inline does improve performance as of Julia 1.10.0 - ((x + 1) / 2)^(2γLoc) * exp(-uLoc * (x+1)^2 / 4) / 2 + ((x + 1) / 2)^(2γ) * exp(-u * (x+1)^2 / 4) / 2 end end f @@ -90,8 +90,8 @@ end function genIntOverlapCore(Δx::T, i₁::Int, α₁::T, i₂::Int, α₂::T, - sameCen::Bool=false) where {T} - getRange = ifelse( sameCen || iszero(Δx), + zeroΔx::Bool=false) where {T} + getRange = ifelse( zeroΔx, Ω::Int -> ifelse(iseven(Ω), begin ΩHalf=Ω÷2; ΩHalf:ΩHalf end, 1:0), Ω::Int -> 0:(Ω÷2) ) f = genIntOverlapTerm(getRange, Δx, i₁, α₁, i₂, α₂) @@ -106,13 +106,13 @@ function ∫overlapCore(::Val{3}, ΔR::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - sameCenAng::NTuple{2, Bool}=(false, false), + sameCen::NTuple{3, Bool}=(false, false, false), coeff::T=T(1.0)) where {T} any(n -> n<0, (ijk₁..., ijk₂...)) && (return T(0.0)) α = α₁ + α₂ res = T(1.0) - for (i₁, i₂, ΔRᵢ) in zip(ijk₁, ijk₂, ΔR) - int = genIntOverlapCore(ΔRᵢ, i₁, α₁, i₂, α₂, sameCenAng[begin]) + for (i₁, i₂, ΔRᵢ, zeroΔx) in zip(ijk₁, ijk₂, ΔR, sameCen) + int = genIntOverlapCore(ΔRᵢ, i₁, α₁, i₂, α₂, zeroΔx) iszero(int) && (return T(0.0)) res *= (-1)^(i₁) * factorial(i₁) * factorial(i₂) * α^(-i₁-i₂) * int end @@ -124,17 +124,16 @@ end R₁::NTuple{3, T}, R₂::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - sameCenAng::NTuple{2, Bool}=(false, false), + sameCen::NTuple{3, Bool}=(false, false, false), coeff::T=T(1.0)) where {T} = -∫overlapCore(Val(3), R₁.-R₂, ijk₁, α₁, ijk₂, α₂, sameCenAng, coeff) +∫overlapCore(Val(3), R₁.-R₂, ijk₁, α₁, ijk₂, α₂, sameCen, coeff) function ∫elecKineticCore(::Val{3}, R₁::NTuple{3, T}, R₂::NTuple{3, T}, - ijk₁::NTuple{3, Int}, α₁::T, + ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - sameCenAng::NTuple{2, Bool}=(false, false)) where {T} - sameCenAng = (sameCenAng[begin], false) + sameCen::NTuple{3, Bool}=(false, false, false)) where {T} ΔR = R₁ .- R₂ shifts = ((1,0,0), (0,1,0), (0,0,1)) mapreduce(+, ijk₁, ijk₂, shifts) do i₁, i₂, Δ𝑙 @@ -142,44 +141,82 @@ function ∫elecKineticCore(::Val{3}, Δijk2 = map(-, ijk₂, Δ𝑙) Δijk3 = map(+, ijk₁, Δ𝑙) Δijk4 = map(+, ijk₂, Δ𝑙) - int1 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk2, α₂, sameCenAng, T(i₁*i₂/2)) - int2 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk4, α₂, sameCenAng, 2α₁*α₂ ) - int3 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk2, α₂, sameCenAng, α₁*i₂ ) - int4 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk4, α₂, sameCenAng, α₂*i₁ ) + int1 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk2, α₂, sameCen, T(i₁*i₂/2)) + int2 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk4, α₂, sameCen, 2α₁*α₂ ) + int3 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk2, α₂, sameCen, α₁*i₂ ) + int4 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk4, α₂, sameCen, α₂*i₁ ) int1 + int2 - int3 - int4 end end function genIntTerm1(Δx::T1, - l₁::T2, o₁::T2, - l₂::T2, o₂::T2, - i₁::T2, α₁::T1, - i₂::T2, α₂::T1) where {T1, T2<:Integer} - (r::T2) -> - ( Δx^muladd(-2, r, o₁+o₂) / (factorial(r) * (factorial∘muladd)(-2, r, o₁+o₂)) ) * - ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * - ( α₂^(o₁-l₂- r) / (factorial(l₂) * factorial(i₂-2l₂-o₂)) ) * - T1( (-1)^(o₂+r) * factorial(o₁+o₂) / - ((1 << 2(l₁+l₂+r)) * factorial(o₁) * factorial(o₂)) ) + l₁::T2, o₁::T2, i₁::T2, α₁::T1, + l₂::T2, o₂::T2, i₂::T2, α₂::T1, zeroΔx::Bool=false) where + {T1, T2<:Integer} + f = let Δx=Δx, l₁=l₁, o₁=o₁, i₁=i₁, α₁=α₁, l₂=l₂, o₂=o₂, i₂=i₂, α₂=α₂, zeroΔx=zeroΔx + if zeroΔx + function (r::T2) + pwr = muladd(-2, r, o₁+o₂) + if iszero(pwr) + ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * + ( α₂^(o₁-l₂- r) / (factorial(l₂) * factorial(i₂-2l₂-o₂)) ) * + ( T1(-1)^(o₂+r) * factorial(o₁+o₂) / ((1 << 2(l₁+l₂+r)) * + factorial(r) * factorial(o₁) * factorial(o₂)) ) + else + T1(0.0) + end + end + else + function (r::T2) + pwr = muladd(-2, r, o₁+o₂) + ( Δx^pwr / (factorial(r ) * factorial(pwr )) ) * + ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * + ( α₂^(o₁-l₂- r) / (factorial(l₂) * factorial(i₂-2l₂-o₂)) ) * + ( T1(-1)^(o₂+r) * factorial(o₁+o₂) / + ((1 << 2(l₁+l₂+r)) * factorial(o₁) * factorial(o₂)) ) + end + end + end + f end -function genIntTerm2core(Δx::T1, μ::T2) where {T1, T2<:Integer} - (u::T2) -> - Δx^(μ-2u) * T1( (-1)^u * factorial(μ) / - ((1 << 2u) * factorial(u) * factorial(μ-2u)) ) +function genIntTerm2core(Δx::T1, μ::T2, zeroΔx::Bool=false) where {T1, T2<:Integer} + f = let Δx=Δx, μ=μ, zeroΔx=zeroΔx + if zeroΔx + function (u::T2) + if μ == 2u + T1(-1)^u * factorial(μ) / ((1 << 2u) * factorial(u)) + else + T1(0.0) + end + end + else + function (u::T2) + Δx^(μ-2u) * ( T1(-1)^u * factorial(μ) / + ((1 << 2u) * factorial(u) * factorial(μ-2u)) ) + end + end + end + f end -function genIntTerm2(Δx::T1, α::T1, o₁::T2, o₂::T2, μ::T2, r::T2) where {T1, T2<:Integer} - (u::T2) -> - genIntTerm2core(Δx, μ)(u) * α^(r-o₁-o₂-u) +function genIntTerm2(Δx::T1, α::T1, o₁::T2, o₂::T2, μ::T2, r::T2, + zeroΔx::Bool=false) where {T1, T2<:Integer} + f = let Δx=Δx, α=α, o₁=o₁, o₂=o₂, μ=μ, r=r, zeroΔx=zeroΔx + function (u::T2) + genIntTerm2core(Δx, μ, zeroΔx)(u) * α^(r-o₁-o₂-u) + end + end + f end function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T, ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - sameCenAng::NTuple{2, Bool}=(false, false)) where {T} + (sameRR₀, sameR₁R₂)::NTuple{2, NTuple{3, Bool}}= + (Tuple∘fill)((false, false, false), 2)) where {T} A = T(0.0) i₁, j₁, k₁ = ijk₁ i₂, j₂, k₂ = ijk₂ @@ -200,13 +237,14 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq₁+opq₂) μsum = sum(μv) Fγs = @inbounds Fγss[μsum+1] - core1s = genIntTerm1.(ΔR₁R₂, lmn₁, opq₁, lmn₂, opq₂, ijk₁, α₁, ijk₂, α₂) + core1s = genIntTerm1.(ΔR₁R₂, lmn₁, opq₁, ijk₁, α₁, lmn₂, opq₂, ijk₂, α₂, + sameR₁R₂) for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) rst = (r, s, t) tmp = T(0.0) - core2s = genIntTerm2.(ΔRR₀, α₁+α₂, opq₁, opq₂, μv, rst) + core2s = genIntTerm2.(ΔRR₀, α₁+α₂, opq₁, opq₂, μv, rst, sameRR₀) for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) γ = μsum - u - v - w @@ -222,8 +260,7 @@ end # function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T, # ijk₁::NTuple{3, Int}, α₁::T, -# ijk₂::NTuple{3, Int}, α₂::T, -# sameCenAng::NTuple{2, Bool}=(false, false)) where {T} +# ijk₂::NTuple{3, Int}, α₂::T) where {T} # A = T(0.0) # i₁, j₁, k₁ = ijk₁ # i₂, j₂, k₂ = ijk₂ @@ -246,7 +283,7 @@ end # μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq₁+opq₂) # μsum = sum(μv) # Fγs = @inbounds Fγss[μsum+1] -# core1s = genIntTerm1.(ΔR₁R₂, lmn₁, opq₁, lmn₂, opq₂, ijk₁, α₁, ijk₂, α₂) +# core1s = genIntTerm1.(ΔR₁R₂, lmn₁, opq₁, ijk₁, α₁, lmn₂, opq₂, ijk₂, α₂) # for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) @@ -274,32 +311,41 @@ end function ∫nucAttractionCore(::Val{3}, Z₀::Int, R₀::NTuple{3, T}, R₁::NTuple{3, T}, R₂::NTuple{3, T}, - ijk₁::NTuple{3, Int}, α₁::T, + ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - sameCenAng::NTuple{2, Bool}=(false, false)) where {T} + sameR₁R₂::NTuple{3, Bool}=(false, false, false)) where {T} α = α₁ + α₂ R = @. (α₁*R₁ + α₂*R₂) / α ΔRR₀ = R .- R₀ ΔR₁R₂ = R₁ .- R₂ β = α * sum(abs2, ΔRR₀) - genIntNucAttCore(ΔRR₀, ΔR₁R₂, β, ijk₁, α₁, ijk₂, α₂, sameCenAng) * + sameRR₀ = iszero.(ΔRR₀) + genIntNucAttCore(ΔRR₀, ΔR₁R₂, β, ijk₁, α₁, ijk₂, α₂, (sameRR₀, sameR₁R₂)) * (π / α) * exp(-α₁ / α * α₂ * sum(abs2, ΔR₁R₂)) * ( -Z₀ * (-1)^sum(ijk₁ .+ ijk₂) * mapMapReduce(ijk₁, factorial) * mapMapReduce(ijk₂, factorial) ) end function genIntTerm3(Δx::T1, - l₁::T2, o₁::T2, - l₂::T2, o₂::T2, - i₁::T2, α₁::T1, - i₂::T2, α₂::T1) where {T1, T2<:Integer} - (r::T2) -> - genIntTerm1(Δx, l₁, o₁, l₂, o₂, i₁, α₁, i₂, α₂)(r) * (α₁+α₂)^muladd(2, l₁+l₂, r) + l₁::T2, o₁::T2, i₁::T2, α₁::T1, + l₂::T2, o₂::T2, i₂::T2, α₂::T1, + zeroΔx::Bool=false) where {T1, T2<:Integer} + f = let Δx=Δx, l₁=l₁, o₁=o₁, i₁=i₁, α₁=α₁, l₂=l₂, o₂=o₂, i₂=i₂, α₂=α₂, zeroΔx=zeroΔx + function (r::T2) + genIntTerm1(Δx, l₁, o₁, i₁, α₁, l₂, o₂, i₂, α₂, zeroΔx)(r) * + (α₁+α₂)^muladd(2, l₁+l₂, r) + end + end + f end -function genIntTerm4(Δx::T1, η::T1, μ::T2) where {T1, T2<:Integer} - (u::T2) -> - genIntTerm2core(Δx, μ)(u) * η^(μ-u) +function genIntTerm4(Δx::T1, η::T1, μ::T2, zeroΔx::Bool=false) where {T1, T2<:Integer} + f = let Δx=Δx, η=η, μ=μ, zeroΔx=zeroΔx + function (u::T2) + genIntTerm2core(Δx, μ, zeroΔx)(u) * η^(μ-u) + end + end + f end @@ -308,7 +354,9 @@ function ∫eeInteractionCore1234(ΔRl::NTuple{3, T}, ΔRr::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, ijk₃::NTuple{3, Int}, α₃::T, - ijk₄::NTuple{3, Int}, α₄::T) where {T} + ijk₄::NTuple{3, Int}, α₄::T, + (sameRl, sameRr, sameRc)::NTuple{3, NTuple{3, Bool}}= + (Tuple∘fill)((false, false, false), 3)) where {T} A = T(0.0) (i₁, j₁, k₁), (i₂, j₂, k₂), (i₃, j₃, k₃), (i₄, j₄, k₄) = ijk₁, ijk₂, ijk₃, ijk₄ @@ -341,9 +389,9 @@ function ∫eeInteractionCore1234(ΔRl::NTuple{3, T}, ΔRr::NTuple{3, T}, μsum = sum(μv) Fγs = F₀toFγ(μsum, β) - core1s = genIntTerm3.(ΔRl, lmn₁, opq₁, lmn₂, opq₂, ijk₁, α₁, ijk₂, α₂) - core2s = genIntTerm3.(ΔRr, lmn₄, opq₄, lmn₃, opq₃, ijk₄, α₄, ijk₃, α₃) - core3s = genIntTerm4.(ΔRc, η, μv) + core1s = genIntTerm3.(ΔRl, lmn₁, opq₁, ijk₁, α₁, lmn₂, opq₂, ijk₂, α₂, sameRl) + core2s = genIntTerm3.(ΔRr, lmn₄, opq₄, ijk₄, α₄, lmn₃, opq₃, ijk₃, α₃, sameRr) + core3s = genIntTerm4.(ΔRc, η, μv, sameRc) for r₁ in 0:((o₁+o₂)÷2), s₁ in 0:((p₁+p₂)÷2), t₁ in 0:((q₁+q₂)÷2), r₂ in 0:((o₃+o₄)÷2), s₂ in 0:((p₃+p₄)÷2), t₂ in 0:((q₃+q₄)÷2) @@ -369,7 +417,9 @@ function ∫eeInteractionCore(::Val{3}, R₁::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, R₂::NTuple{3, T}, ijk₂::NTuple{3, Int}, α₂::T, R₃::NTuple{3, T}, ijk₃::NTuple{3, Int}, α₃::T, - R₄::NTuple{3, T}, ijk₄::NTuple{3, Int}, α₄::T) where {T} + R₄::NTuple{3, T}, ijk₄::NTuple{3, Int}, α₄::T, + (sameR₁R₂, sameR₃R₄)::NTuple{2, NTuple{3, Bool}}= + (Tuple∘fill)((false, false, false), 2)) where {T} ΔRl = R₁ .- R₂ ΔRr = R₃ .- R₄ αl = α₁ + α₂ @@ -377,9 +427,11 @@ function ∫eeInteractionCore(::Val{3}, ηl = α₁ / αl * α₂ ηr = α₃ / αr * α₄ ΔRc = @. (α₁*R₁ + α₂*R₂) / αl - (α₃*R₃ + α₄*R₄) / αr + zeroΔRc = iszero.(ΔRc) η = αl / (α₁ + α₂ + α₃ + α₄) * αr β = η * sum(abs2, ΔRc) - ∫eeInteractionCore1234(ΔRl, ΔRr, ΔRc, β, η, ijk₁, α₁, ijk₂, α₂, ijk₃, α₃, ijk₄, α₄) * + ∫eeInteractionCore1234(ΔRl, ΔRr, ΔRc, β, η, ijk₁, α₁, ijk₂, α₂, ijk₃, α₃, ijk₄, α₄, + (sameR₁R₂, sameR₃R₄, zeroΔRc)) * T(πvals[2.5]) / (αl * αr * sqrt(αl + αr)) * exp(-ηl * sum(abs2, ΔRl)) * exp(-ηr * sum(abs2, ΔRr)) * mapreduce(*, ijk₁, ijk₂, ijk₃, ijk₄) do l₁, l₂, l₃, l₄ @@ -624,10 +676,10 @@ function getOneBodyInt(::Type{T}, ::Val{D}, ∫1e::F, @nospecialize(optPosArgs:: {T, D, F<:Function} (R₁, ijk₁, ps₁, 𝑙₁), (R₂, ijk₂, ps₂, 𝑙₂) = reformatIntData1(iBl, bfs) 𝑙₁==𝑙₂==0 || isIntZero(F, optPosArgs, R₁,R₂, ijk₁,ijk₂) && (return T(0.0)) - sameCenAng = (R₁==R₂, ijk₁==ijk₂) - uniquePairs, uPairCoeffs = getOneBodyUniquePairs(prod(sameCenAng), ps₁, ps₂) + sameCen = iszero.(R₁ .- R₂) + uniquePairs, uPairCoeffs = getOneBodyUniquePairs(prod(sameCen)*(ijk₁==ijk₂), ps₁, ps₂) mapreduce(+, uniquePairs, uPairCoeffs) do x, y - ∫1e(Val(D), optPosArgs..., R₁, R₂, ijk₁, x[1], ijk₂, x[2], sameCenAng)::T * y + ∫1e(Val(D), optPosArgs..., R₁, R₂, ijk₁, x[1], ijk₂, x[2], sameCen)::T * y end end @@ -988,8 +1040,12 @@ function getTwoBodyInt(::Type{T}, ::Val{D}, ∫2e::F, @nospecialize(optPosArgs:: 𝑙₁==𝑙₂==𝑙₃==𝑙₄==0 || isIntZero(F, optPosArgs, R₁,R₂,R₃,R₄, ijk₁,ijk₂,ijk₃,ijk₄) && (return T(0.0)) - f1 = (R₁ == R₂ && ijk₁ == ijk₂) - f2 = (R₃ == R₄ && ijk₃ == ijk₄) + sameR₁R₂ = iszero.(R₁ .- R₂) + sameR₃R₄ = iszero.(R₃ .- R₄) + sameCens = (sameR₁R₂, sameR₃R₄) + + f1 = (prod(sameR₁R₂) && ijk₁ == ijk₂) + f2 = (prod(sameR₃R₄) && ijk₃ == ijk₄) f3 = (R₁ == R₃ && ijk₁ == ijk₃ && R₂ == R₄ && ijk₂ == ijk₄) f4 = (R₁ == R₄ && ijk₁ == ijk₄) f5 = (R₂ == R₃ && ijk₂ == ijk₃) @@ -1000,7 +1056,7 @@ function getTwoBodyInt(::Type{T}, ::Val{D}, ∫2e::F, @nospecialize(optPosArgs:: ∫2e(Val(D), optPosArgs..., R₁, ijk₁, x[1], R₂, ijk₂, x[2], R₃, ijk₃, x[3], - R₄, ijk₄, x[4])::T * y + R₄, ijk₄, x[4], sameCens)::T * y end |> sum # Fewer allocations than mapreduce. end From 436809b94b296ff57b4e5d57e4701448dc6431ad Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Wed, 13 Mar 2024 17:44:02 -0400 Subject: [PATCH 12/22] Code optimization. --- src/Tools.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Tools.jl b/src/Tools.jl index 4b19ba19..ee7363d2 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -994,9 +994,9 @@ end mapMapReduce(tp::NTuple{N}, fs::NTuple{N, Function}, op::F=*) where {N, F} = -mapreduce((x, y)->y(x), op, tp, fs) +reduce(op, map((x, y)->y(x), tp, fs)) -mapMapReduce(tp::NTuple{N}, f::F1, op::F2=*) where {N, F1, F2} = mapreduce(f, op, tp) +mapMapReduce(tp::NTuple{N}, f::F1, op::F2=*) where {N, F1, F2} = reduce(op, map(f, tp)) rmsOf(arr::AbstractArray) = norm(arr) / (sqrt∘length)(arr) From e3b003de8681f435a3bb2edeb04860ef82ed3301 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Thu, 14 Mar 2024 01:29:42 -0400 Subject: [PATCH 13/22] Code optimization. --- src/Integrals/Core.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index ff4b2e1d..008a45fe 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -26,7 +26,7 @@ const MaxGQpointNum = 10000 function FγCore(γ::Int, u::T, GQpointNum::Int) where {T} GQpointNum = min(GQpointNum, MaxGQpointNum) - GQnodes, GQweights = get!(NodesAndWeightsOfGQ, GQpointNum) do + GQnodes, GQweights = LRUCache.get!(NodesAndWeightsOfGQ, GQpointNum) do gausslegendre(GQpointNum) end GQnodes = convert(Vector{T}, GQnodes) @@ -176,6 +176,8 @@ function genIntTerm1(Δx::T1, ( T1(-1)^(o₂+r) * factorial(o₁+o₂) / ((1 << 2(l₁+l₂+r)) * factorial(o₁) * factorial(o₂)) ) end + + # Take T1(-1)^(o₂+r) outside of the closure to make it symmetric end end f @@ -678,9 +680,9 @@ function getOneBodyInt(::Type{T}, ::Val{D}, ∫1e::F, @nospecialize(optPosArgs:: 𝑙₁==𝑙₂==0 || isIntZero(F, optPosArgs, R₁,R₂, ijk₁,ijk₂) && (return T(0.0)) sameCen = iszero.(R₁ .- R₂) uniquePairs, uPairCoeffs = getOneBodyUniquePairs(prod(sameCen)*(ijk₁==ijk₂), ps₁, ps₂) - mapreduce(+, uniquePairs, uPairCoeffs) do x, y + map(uniquePairs, uPairCoeffs) do x, y ∫1e(Val(D), optPosArgs..., R₁, R₂, ijk₁, x[1], ijk₂, x[2], sameCen)::T * y - end + end |> sum end @@ -1093,10 +1095,10 @@ get1BCompInt(::Type{T}, ::Val{D}, ::typeof(∫nucAttractionCore), Tuple{NTuple{D, T}, Vararg{NTuple{D, T}, NNMO}}}, iBl::Union{iBlTs[1], iBlTs[3]}, ::NTuple{2, Int}, bfs::NTupleOfFGTBF{2, T, D}) where {T, D, NNMO} = -mapreduce(+, nucAndCoords[1], nucAndCoords[2]) do ele, coord +map(nucAndCoords[1], nucAndCoords[2]) do ele, coord getOneBodyInt(T, Val(D), ∫nucAttractionCore, (getCharge(ele), coord), iBl, orderFGTBG(bfs)) -end +end |> sum # j==i j!=i const Int1eBIndexLabels = Dict([( true,), (false,)] .=> [Val(:aa), Val(:ab)]) From d2a7c2648671ae7f53dc08993541975d4e15cdf2 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Sat, 16 Mar 2024 16:41:31 -0400 Subject: [PATCH 14/22] Code optimization. --- src/Integrals/Core.jl | 47 ++++++++++++++++++------------------------- src/Tools.jl | 7 +++++-- 2 files changed, 25 insertions(+), 29 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 008a45fe..3d250255 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -78,9 +78,8 @@ function genIntOverlapTerm(getRangeFunc::F, Δx::T, i₁::Int, α₁::T, res += Δx^(Ω-2o) * ( α₁^(i₂ - l₁ - 2l₂ - o) / (factorial(l₂) * factorial(i₁-2l₁)) ) * ( α₂^(i₁ - l₂ - 2l₁ - o) / (factorial(l₁) * factorial(i₂-2l₂)) ) * - ( T(-1)^o * factorial(Ω) / - ((1 << 2(l₁ + l₂ + o)) * factorial(o ) * factorial(Ω -2o )) ) * - (α₁ + α₂)^muladd(2, (l₁ + l₂), o) + ( n1Power(o) * factorial(Ω) * (α₁ + α₂)^muladd(2, (l₁ + l₂), o) / + ((1 << 2(l₁ + l₂ + o)) * factorial(o ) * factorial(Ω -2o )) ) end res end @@ -114,7 +113,7 @@ function ∫overlapCore(::Val{3}, for (i₁, i₂, ΔRᵢ, zeroΔx) in zip(ijk₁, ijk₂, ΔR, sameCen) int = genIntOverlapCore(ΔRᵢ, i₁, α₁, i₂, α₂, zeroΔx) iszero(int) && (return T(0.0)) - res *= (-1)^(i₁) * factorial(i₁) * factorial(i₂) * α^(-i₁-i₂) * int + res *= n1Power(i₁) * factorial(i₁) * factorial(i₂) * α^(-i₁-i₂) * int end res *= sqrt((π/α)^3) * exp(-α₁ / α * α₂ * sum(abs2, ΔR)) res * coeff @@ -150,34 +149,27 @@ function ∫elecKineticCore(::Val{3}, end +function genIntTerm1Core(r::Int, l₁::Int, o₁::Int, i₁::Int, α₁::T, + l₂::Int, o₂::Int, i₂::Int, α₂::T) where {T} + ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * + ( α₂^(o₁-l₂- r) / (factorial(l₂) * factorial(i₂-2l₂-o₂)) ) * + ( n1Power(o₂+r, T) * factorial(o₁+o₂) / ((1 << 2(l₁+l₂+r)) * + factorial(r) * factorial(o₁) * factorial(o₂)) ) +end + function genIntTerm1(Δx::T1, l₁::T2, o₁::T2, i₁::T2, α₁::T1, l₂::T2, o₂::T2, i₂::T2, α₂::T1, zeroΔx::Bool=false) where {T1, T2<:Integer} f = let Δx=Δx, l₁=l₁, o₁=o₁, i₁=i₁, α₁=α₁, l₂=l₂, o₂=o₂, i₂=i₂, α₂=α₂, zeroΔx=zeroΔx - if zeroΔx - function (r::T2) - pwr = muladd(-2, r, o₁+o₂) - if iszero(pwr) - ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * - ( α₂^(o₁-l₂- r) / (factorial(l₂) * factorial(i₂-2l₂-o₂)) ) * - ( T1(-1)^(o₂+r) * factorial(o₁+o₂) / ((1 << 2(l₁+l₂+r)) * - factorial(r) * factorial(o₁) * factorial(o₂)) ) - else - T1(0.0) - end - end - else - function (r::T2) - pwr = muladd(-2, r, o₁+o₂) - ( Δx^pwr / (factorial(r ) * factorial(pwr )) ) * - ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * - ( α₂^(o₁-l₂- r) / (factorial(l₂) * factorial(i₂-2l₂-o₂)) ) * - ( T1(-1)^(o₂+r) * factorial(o₁+o₂) / - ((1 << 2(l₁+l₂+r)) * factorial(o₁) * factorial(o₂)) ) + function (r::T2) + pwr = muladd(-2, r, o₁+o₂) + coeff = if zeroΔx + iszero(pwr) ? T1(1) : (return T1(0.0)) + else + Δx^pwr / factorial(pwr) end - - # Take T1(-1)^(o₂+r) outside of the closure to make it symmetric + coeff * genIntTerm1Core(r, l₁, o₁, i₁, α₁, l₂, o₂, i₂, α₂) end end f @@ -249,7 +241,8 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T core2s = genIntTerm2.(ΔRR₀, α₁+α₂, opq₁, opq₂, μv, rst, sameRR₀) for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) - γ = μsum - u - v - w + uvw = u + v + w + γ = μsum - uvw @inbounds tmp += mapMapReduce((u, v, w), core2s) * 2Fγs[γ+1] end A += mapMapReduce(rst, core1s) * tmp diff --git a/src/Tools.jl b/src/Tools.jl index ee7363d2..7df942ed 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -878,7 +878,7 @@ function skipIndices(arr::AbstractArray{Int}, ints::AbstractVector{Int}) else all(i > 0 for i in ints) || throw(DomainError(ints, "Every element of `ints` should be positive.")) - maxIdx = max(arr...) + maxIdx = maximum(arr) maxIdxN = maxIdx + length(ints) ints = filter!(x->x<=maxIdxN, sort(ints)) idsN = deleteat!(collect(1:maxIdxN), ints) @@ -1027,4 +1027,7 @@ function betterSum(x) # Improved Kahan–Babuška algorithm fro array summation s = t end s .+ c -end \ No newline at end of file +end + +n1Power(x::Int) = ifelse(isodd(x), -1, 1) +n1Power(x::Int, ::Type{T}) where {T} = ifelse(isodd(x), T(-1), T(1)) \ No newline at end of file From 5adf6c84c9e68eb0c6ee64bd65076a4c3eb3f4d6 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Sat, 16 Mar 2024 16:42:11 -0400 Subject: [PATCH 15/22] Replaced `max` with `maximum`. --- src/HartreeFock.jl | 2 +- test/test-functions/Shared.jl | 2 +- test/unit-tests/Basis-test.jl | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/HartreeFock.jl b/src/HartreeFock.jl index 6b54839b..7398af07 100644 --- a/src/HartreeFock.jl +++ b/src/HartreeFock.jl @@ -888,7 +888,7 @@ function runHFcore(bs::GTBasis{T, D, BN, BFT}, Ntot = (N isa Int) ? N : (N[begin] + N[end]) Ntot > Nlow || throw(DomainError(N, "$(HFT) requires more than $(Nlow) electrons.")) Ns = splitSpins(Val(HFT), N) - leastNb = max(Ns...) + leastNb = maximum(Ns) BN < leastNb && throw(DomainError(BN, "The number of basis functions should be no "* "less than $(leastNb).")) nuc = arrayToTuple(nuc) diff --git a/test/test-functions/Shared.jl b/test/test-functions/Shared.jl index d05bb2fd..313cb074 100644 --- a/test/test-functions/Shared.jl +++ b/test/test-functions/Shared.jl @@ -58,7 +58,7 @@ function compr2Arrays3(cprTuple::NamedTuple{<:Any, <:Tuple{T1, T2}}, atol::Real, showAllDiff && !res && println("$(name1) - $(name2) = ", diff) println(additionalInfo) v, i = findmax(abs, diff) - println("max(abs.($(name1) - $(name2))...) = ", v, " index = ", i) + println("maximum(abs.($(name1) - $(name2))) = ", v, " index = ", i) end res end \ No newline at end of file diff --git a/test/unit-tests/Basis-test.jl b/test/unit-tests/Basis-test.jl index 596c174e..e6f56892 100644 --- a/test/unit-tests/Basis-test.jl +++ b/test/unit-tests/Basis-test.jl @@ -957,9 +957,9 @@ bfmN = normalizeBasis(bfm) @test isapprox(overlap(bfmN, bfmN), 1, atol=t) bf4_3N = normalizeBasis(bf4_3) @test typeof(bf4_3N[]) == typeof(bf4_3) -@test isapprox(max( abs.((overlaps(bf4_3N)|>diag) .- 1)... ), 0, atol=t) +@test isapprox(maximum( abs.((overlaps(bf4_3N)|>diag) .- 1) ), 0, atol=t) bf4_4N = normalizeBasis(bf4_4) @test length(bf4_4N) == 2 -@test isapprox(max( abs.((overlaps(bf4_4N)|>diag) .- 1)... ), 0, atol=t) +@test isapprox(maximum( abs.((overlaps(bf4_4N)|>diag) .- 1) ), 0, atol=t) end \ No newline at end of file From accd993fe4e16f42231cd4a18128c3b63505df1c Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Sun, 17 Mar 2024 03:57:06 -0400 Subject: [PATCH 16/22] Optimized integral-term computation redundancy. --- src/Integrals/Core.jl | 228 +++++++++++++++++------------------------- 1 file changed, 93 insertions(+), 135 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 3d250255..9d558ec3 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -149,60 +149,44 @@ function ∫elecKineticCore(::Val{3}, end -function genIntTerm1Core(r::Int, l₁::Int, o₁::Int, i₁::Int, α₁::T, - l₂::Int, o₂::Int, i₂::Int, α₂::T) where {T} - ( α₁^(o₂-l₁- r) / (factorial(l₁) * factorial(i₁-2l₁-o₁)) ) * - ( α₂^(o₁-l₂- r) / (factorial(l₂) * factorial(i₂-2l₂-o₂)) ) * - ( n1Power(o₂+r, T) * factorial(o₁+o₂) / ((1 << 2(l₁+l₂+r)) * - factorial(r) * factorial(o₁) * factorial(o₂)) ) -end - -function genIntTerm1(Δx::T1, - l₁::T2, o₁::T2, i₁::T2, α₁::T1, - l₂::T2, o₂::T2, i₂::T2, α₂::T1, zeroΔx::Bool=false) where - {T1, T2<:Integer} - f = let Δx=Δx, l₁=l₁, o₁=o₁, i₁=i₁, α₁=α₁, l₂=l₂, o₂=o₂, i₂=i₂, α₂=α₂, zeroΔx=zeroΔx - function (r::T2) - pwr = muladd(-2, r, o₁+o₂) +function genIntTerm(arg1::T1, Δx::T2, zeroΔx::Bool) where {T1<:Integer, T2<:Real} + fRes = let arg1=arg1, Δx=Δx, zeroΔx=zeroΔx + function (arg2::T1) + pwr = muladd(-2, arg2, arg1) coeff = if zeroΔx - iszero(pwr) ? T1(1) : (return T1(0.0)) + iszero(pwr) ? T2(1) : (return T2(0)) else Δx^pwr / factorial(pwr) end - coeff * genIntTerm1Core(r, l₁, o₁, i₁, α₁, l₂, o₂, i₂, α₂) + coeff / factorial(arg2) end end - f + fRes end -function genIntTerm2core(Δx::T1, μ::T2, zeroΔx::Bool=false) where {T1, T2<:Integer} - f = let Δx=Δx, μ=μ, zeroΔx=zeroΔx - if zeroΔx - function (u::T2) - if μ == 2u - T1(-1)^u * factorial(μ) / ((1 << 2u) * factorial(u)) - else - T1(0.0) - end - end - else - function (u::T2) - Δx^(μ-2u) * ( T1(-1)^u * factorial(μ) / - ((1 << 2u) * factorial(u) * factorial(μ-2u)) ) - end - end - end - f +function termProd1(lmn₁::NTuple{3, T1}, lmn₂::NTuple{3, T1}) where {T1<:Integer} + prod( @. factorial(lmn₁) * factorial(lmn₂) ) end -function genIntTerm2(Δx::T1, α::T1, o₁::T2, o₂::T2, μ::T2, r::T2, - zeroΔx::Bool=false) where {T1, T2<:Integer} - f = let Δx=Δx, α=α, o₁=o₁, o₂=o₂, μ=μ, r=r, zeroΔx=zeroΔx - function (u::T2) - genIntTerm2core(Δx, μ, zeroΔx)(u) * α^(r-o₁-o₂-u) - end - end - f +function termProd2(v::NTuple{3, T1}) where {T1<:Integer} + prod( factorial.(v) ) +end + +function termProd3(ijk₁::NTuple{3, T1}, lmn₁::NTuple{3, T1}, opq₁::NTuple{3, T1}, + ijk₂::NTuple{3, T1}, lmn₂::NTuple{3, T1}, opq₂::NTuple{3, T1}) where + {T1<:Integer} + prod( @. ( factorial(opq₁) * factorial(ijk₁-2lmn₁-opq₁) * + factorial(opq₂) * factorial(ijk₂-2lmn₂-opq₂) ) ) +end +function termProd4(lmn₁Sum::T1, opq₁Sum::T1, lmn₂Sum::T1, opq₂Sum::T1, rstSum::T1, + (α₁, α₂)::NTuple{2, T2}) where {T1<:Integer, T2<:Real} + ( n1Power(opq₂Sum + rstSum) * α₁^(opq₂Sum - lmn₁Sum - rstSum) * + α₂^(opq₁Sum - lmn₂Sum - rstSum) ) / + ( 1<<2(lmn₁Sum + lmn₂Sum + rstSum) ) +end + +function termProd5(uvwSum::T1, ::Type{T2}) where {T1<:Integer, T2<:Real} + n1Power(uvwSum, T2) / (1 << 2uvwSum) end @@ -212,6 +196,7 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T (sameRR₀, sameR₁R₂)::NTuple{2, NTuple{3, Bool}}= (Tuple∘fill)((false, false, false), 2)) where {T} A = T(0.0) + α = α₁ + α₂ i₁, j₁, k₁ = ijk₁ i₂, j₂, k₂ = ijk₂ ijkSum = sum(ijk₁) + sum(ijk₂) @@ -221,31 +206,46 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T lmn₁ = (l₁, m₁, n₁) lmn₂ = (l₂, m₂, n₂) + lmn₁Sum = sum(lmn₁) + lmn₂Sum = sum(lmn₂) + term1Coeff1 = termProd1(lmn₁, lmn₂) for o₁ in 0:(i₁-2l₁), p₁ in 0:(j₁-2m₁), q₁ in 0:(k₁-2n₁), o₂ in 0:(i₂-2l₂), p₂ in 0:(j₂-2m₂), q₂ in 0:(k₂-2n₂) opq₁ = (o₁, p₁, q₁) opq₂ = (o₂, p₂, q₂) - - μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq₁+opq₂) - μsum = sum(μv) - Fγs = @inbounds Fγss[μsum+1] - core1s = genIntTerm1.(ΔR₁R₂, lmn₁, opq₁, ijk₁, α₁, lmn₂, opq₂, ijk₂, α₂, - sameR₁R₂) + opq = opq₁ .+ opq₂ + opq₁Sum = sum(opq₁) + opq₂Sum = sum(opq₂) + opqSum = opq₁Sum + opq₂Sum + term1Coeff2 = termProd2(opq) + term1Coeff3 = termProd3(ijk₁, lmn₁, opq₁, ijk₂, lmn₂, opq₂) + + μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq) + term2Coeff1 = termProd2(μv) + μSum = sum(μv) + Fγs = @inbounds Fγss[μSum+1] + core1s = genIntTerm.(opq₁.+opq₂, ΔR₁R₂, sameR₁R₂) for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) rst = (r, s, t) + rstSum = sum(rst) + term1Coeff4 = termProd4(lmn₁Sum, opq₁Sum, lmn₂Sum, opq₂Sum, rstSum, (α₁,α₂)) tmp = T(0.0) - core2s = genIntTerm2.(ΔRR₀, α₁+α₂, opq₁, opq₂, μv, rst, sameRR₀) + core2s = genIntTerm.(μv, ΔRR₀, sameRR₀) for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) - uvw = u + v + w - γ = μsum - uvw - @inbounds tmp += mapMapReduce((u, v, w), core2s) * 2Fγs[γ+1] + uvwSum = u + v + w + γ = μSum - uvwSum + term2Coeff2 = termProd5(uvwSum, T) * 2Fγs[γ+1] * + α^(rstSum - opqSum - uvwSum) + @inbounds tmp += mapMapReduce((u, v, w), core2s) * + term2Coeff1 * term2Coeff2 end - A += mapMapReduce(rst, core1s) * tmp + A += mapMapReduce(rst, core1s) * tmp * + (term1Coeff2 * term1Coeff4 / (term1Coeff1 * term1Coeff3)) end end @@ -253,56 +253,6 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T A end -# function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T, -# ijk₁::NTuple{3, Int}, α₁::T, -# ijk₂::NTuple{3, Int}, α₂::T) where {T} -# A = T(0.0) -# i₁, j₁, k₁ = ijk₁ -# i₂, j₂, k₂ = ijk₂ -# ijkSum = sum(ijk₁) + sum(ijk₂) -# Fγss = [F₀toFγ(γ, β) for γ in 0:ijkSum] -# for l₁ in 0:(i₁÷2), m₁ in 0:(j₁÷2), n₁ in 0:(k₁÷2), -# l₂ in 0:(i₂÷2), m₂ in 0:(j₂÷2), n₂ in 0:(k₂÷2) - -# lmn₁ = (l₁, m₁, n₁) -# lmn₂ = (l₂, m₂, n₂) - -# res = T(0) - -# for o₁ in 0:(i₁-2l₁), p₁ in 0:(j₁-2m₁), q₁ in 0:(k₁-2n₁), -# o₂ in 0:(i₂-2l₂), p₂ in 0:(j₂-2m₂), q₂ in 0:(k₂-2n₂) - -# opq₁ = (o₁, p₁, q₁) -# opq₂ = (o₂, p₂, q₂) - -# μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq₁+opq₂) -# μsum = sum(μv) -# Fγs = @inbounds Fγss[μsum+1] -# core1s = genIntTerm1.(ΔR₁R₂, lmn₁, opq₁, ijk₁, α₁, lmn₂, opq₂, ijk₂, α₂) - -# for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) - -# rst = (r, s, t) -# tmp = T(0.0) -# core2s = genIntTerm2.(ΔRR₀, α₁+α₂, opq₁, opq₂, μv, rst) - -# for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) -# γ = μsum - u - v - w -# @inbounds tmp += mapMapReduce((u, v, w), core2s) * 2Fγs[γ+1] -# end -# rr1 = mapMapReduce(rst, core1s) -# rr2 = tmp -# println(((rr1, (o₁+o₂)÷2, (p₁+p₂)÷2, (q₁+q₂)÷2), (rr2, μˣ, μʸ, μᶻ)), ",") -# res += rr1 * rr2 -# end -# end -# A += res -# end -# println() -# A -# end - - function ∫nucAttractionCore(::Val{3}, Z₀::Int, R₀::NTuple{3, T}, R₁::NTuple{3, T}, R₂::NTuple{3, T}, @@ -321,28 +271,6 @@ function ∫nucAttractionCore(::Val{3}, mapMapReduce(ijk₁, factorial) * mapMapReduce(ijk₂, factorial) ) end -function genIntTerm3(Δx::T1, - l₁::T2, o₁::T2, i₁::T2, α₁::T1, - l₂::T2, o₂::T2, i₂::T2, α₂::T1, - zeroΔx::Bool=false) where {T1, T2<:Integer} - f = let Δx=Δx, l₁=l₁, o₁=o₁, i₁=i₁, α₁=α₁, l₂=l₂, o₂=o₂, i₂=i₂, α₂=α₂, zeroΔx=zeroΔx - function (r::T2) - genIntTerm1(Δx, l₁, o₁, i₁, α₁, l₂, o₂, i₂, α₂, zeroΔx)(r) * - (α₁+α₂)^muladd(2, l₁+l₂, r) - end - end - f -end - -function genIntTerm4(Δx::T1, η::T1, μ::T2, zeroΔx::Bool=false) where {T1, T2<:Integer} - f = let Δx=Δx, η=η, μ=μ, zeroΔx=zeroΔx - function (u::T2) - genIntTerm2core(Δx, μ, zeroΔx)(u) * η^(μ-u) - end - end - f -end - function ∫eeInteractionCore1234(ΔRl::NTuple{3, T}, ΔRr::NTuple{3, T}, ΔRc::NTuple{3, T}, β::T, η::T, @@ -356,6 +284,8 @@ function ∫eeInteractionCore1234(ΔRl::NTuple{3, T}, ΔRr::NTuple{3, T}, (i₁, j₁, k₁), (i₂, j₂, k₂), (i₃, j₃, k₃), (i₄, j₄, k₄) = ijk₁, ijk₂, ijk₃, ijk₄ IJK = @. ijk₁ + ijk₂ + ijk₃ + ijk₄ + αL = α₁ + α₂ + αR = α₃ + α₄ for l₁ in 0:(i₁÷2), m₁ in 0:(j₁÷2), n₁ in 0:(k₁÷2), l₂ in 0:(i₂÷2), m₂ in 0:(j₂÷2), n₂ in 0:(k₂÷2), @@ -366,6 +296,14 @@ function ∫eeInteractionCore1234(ΔRl::NTuple{3, T}, ΔRr::NTuple{3, T}, lmn₂ = (l₂, m₂, n₂) lmn₃ = (l₃, m₃, n₃) lmn₄ = (l₄, m₄, n₄) + lmn₁Sum = sum(lmn₁) + lmn₂Sum = sum(lmn₂) + lmn₃Sum = sum(lmn₃) + lmn₄Sum = sum(lmn₄) + lmnLSum = lmn₁Sum + lmn₂Sum + lmnRSum = lmn₃Sum + lmn₄Sum + term1Coeff1 = termProd1(lmn₁, lmn₂) + term2Coeff1 = termProd1(lmn₄, lmn₃) for o₁ in 0:(i₁-2l₁), p₁ in 0:(j₁-2m₁), q₁ in 0:(k₁-2n₁), o₂ in 0:(i₂-2l₂), p₂ in 0:(j₂-2m₂), q₂ in 0:(k₂-2n₂), @@ -376,30 +314,50 @@ function ∫eeInteractionCore1234(ΔRl::NTuple{3, T}, ΔRr::NTuple{3, T}, opq₂ = (o₂, p₂, q₂) opq₃ = (o₃, p₃, q₃) opq₄ = (o₄, p₄, q₄) + opqL = opq₁ .+ opq₂ + opqR = opq₃ .+ opq₄ + opq₁Sum = sum(opq₁) + opq₂Sum = sum(opq₂) + opq₃Sum = sum(opq₃) + opq₄Sum = sum(opq₄) + term1Coeff2 = termProd2(opqL) + term2Coeff2 = termProd2(opqR) + term1Coeff3 = termProd3(ijk₁, lmn₁, opq₁, ijk₂, lmn₂, opq₂) + term2Coeff3 = termProd3(ijk₄, lmn₄, opq₄, ijk₃, lmn₃, opq₃) μˣ, μʸ, μᶻ = μv = begin - @. IJK - muladd(2, lmn₁+lmn₂+lmn₃+lmn₄, opq₁+opq₂+opq₃+opq₄) + @. IJK - muladd(2, lmn₁+lmn₂+lmn₃+lmn₄, opqL+opqR) end - + term3Coeff1 = termProd2(μv) μsum = sum(μv) Fγs = F₀toFγ(μsum, β) - core1s = genIntTerm3.(ΔRl, lmn₁, opq₁, ijk₁, α₁, lmn₂, opq₂, ijk₂, α₂, sameRl) - core2s = genIntTerm3.(ΔRr, lmn₄, opq₄, ijk₄, α₄, lmn₃, opq₃, ijk₃, α₃, sameRr) - core3s = genIntTerm4.(ΔRc, η, μv, sameRc) + core1s = genIntTerm.(opqL, ΔRl, sameRl) + core2s = genIntTerm.(opqR, ΔRr, sameRr) + core3s = genIntTerm.(μv, ΔRc, sameRc) for r₁ in 0:((o₁+o₂)÷2), s₁ in 0:((p₁+p₂)÷2), t₁ in 0:((q₁+q₂)÷2), r₂ in 0:((o₃+o₄)÷2), s₂ in 0:((p₃+p₄)÷2), t₂ in 0:((q₃+q₄)÷2) rst₁ = (r₁, s₁, t₁) rst₂ = (r₂, s₂, t₂) + rst₁Sum = sum(rst₁) + rst₂Sum = sum(rst₂) + term1Coeff4 = termProd4(lmn₁Sum, opq₁Sum, lmn₂Sum, opq₂Sum, rst₁Sum, (α₁, α₂)) * (αL)^muladd(2, lmnLSum, rst₁Sum) + term2Coeff4 = termProd4(lmn₄Sum, opq₄Sum, lmn₃Sum, opq₃Sum, rst₂Sum, (α₄, α₃)) * (αR)^muladd(2, lmnRSum, rst₂Sum) tmp = T(0.0) for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) - γ = μsum - u - v - w - @inbounds tmp += mapMapReduce((u, v, w), core3s) * 2Fγs[γ+1] + uvwSum = u + v + w + γ = μsum - uvwSum + term3Coeff2 = termProd5(uvwSum, T) * 2Fγs[γ+1] * η^γ + @inbounds tmp += mapMapReduce((u, v, w), core3s) * + term3Coeff1 * term3Coeff2 end - A += mapMapReduce(rst₁, core1s) * mapMapReduce(rst₂, core2s) * tmp + A += mapMapReduce(rst₁, core1s) * mapMapReduce(rst₂, core2s) * + (term1Coeff2 * term1Coeff4 / (term1Coeff1 * term1Coeff3)) * + (term2Coeff2 * term2Coeff4 / (term2Coeff1 * term2Coeff3)) * + tmp end end From c0e483609e0f84eb940834d6377841a7ae00857b Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Sun, 17 Mar 2024 16:51:47 -0400 Subject: [PATCH 17/22] Code formatting. --- src/Integrals/Core.jl | 73 +++++++++++++++++++++---------------------- 1 file changed, 36 insertions(+), 37 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 9d558ec3..c7676fa7 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -154,7 +154,7 @@ function genIntTerm(arg1::T1, Δx::T2, zeroΔx::Bool) where {T1<:Integer, T2<:Re function (arg2::T1) pwr = muladd(-2, arg2, arg1) coeff = if zeroΔx - iszero(pwr) ? T2(1) : (return T2(0)) + iszero(pwr) ? T2(1.0) : (return T2(0.0)) else Δx^pwr / factorial(pwr) end @@ -173,8 +173,9 @@ function termProd2(v::NTuple{3, T1}) where {T1<:Integer} end function termProd3(ijk₁::NTuple{3, T1}, lmn₁::NTuple{3, T1}, opq₁::NTuple{3, T1}, - ijk₂::NTuple{3, T1}, lmn₂::NTuple{3, T1}, opq₂::NTuple{3, T1}) where - {T1<:Integer} + ijk₂::NTuple{3, T1}, lmn₂::NTuple{3, T1}, opq₂::NTuple{3, T1}, + ::Type{T2}) where {T1<:Integer, T2<:Real} + (T2∘termProd2)(opq₁ .+ opq₂) / prod( @. ( factorial(opq₁) * factorial(ijk₁-2lmn₁-opq₁) * factorial(opq₂) * factorial(ijk₂-2lmn₂-opq₂) ) ) end @@ -208,7 +209,7 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T lmn₂ = (l₂, m₂, n₂) lmn₁Sum = sum(lmn₁) lmn₂Sum = sum(lmn₂) - term1Coeff1 = termProd1(lmn₁, lmn₂) + subterm1 = termProd1(lmn₁, lmn₂) for o₁ in 0:(i₁-2l₁), p₁ in 0:(j₁-2m₁), q₁ in 0:(k₁-2n₁), o₂ in 0:(i₂-2l₂), p₂ in 0:(j₂-2m₂), q₂ in 0:(k₂-2n₂) @@ -219,33 +220,31 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T opq₁Sum = sum(opq₁) opq₂Sum = sum(opq₂) opqSum = opq₁Sum + opq₂Sum - term1Coeff2 = termProd2(opq) - term1Coeff3 = termProd3(ijk₁, lmn₁, opq₁, ijk₂, lmn₂, opq₂) + subterm2 = termProd3(ijk₁, lmn₁, opq₁, ijk₂, lmn₂, opq₂, T) / subterm1 μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq) - term2Coeff1 = termProd2(μv) + subterm3 = termProd2(μv) μSum = sum(μv) Fγs = @inbounds Fγss[μSum+1] - core1s = genIntTerm.(opq₁.+opq₂, ΔR₁R₂, sameR₁R₂) + genTerm1s = genIntTerm.(opq₁.+opq₂, ΔR₁R₂, sameR₁R₂) for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) rst = (r, s, t) rstSum = sum(rst) - term1Coeff4 = termProd4(lmn₁Sum, opq₁Sum, lmn₂Sum, opq₂Sum, rstSum, (α₁,α₂)) - tmp = T(0.0) - core2s = genIntTerm.(μv, ΔRR₀, sameRR₀) + term1Coeff = termProd4(lmn₁Sum, opq₁Sum, lmn₂Sum, opq₂Sum, rstSum, + (α₁, α₂)) * subterm2 + term2 = T(0.0) + genTerm2s = genIntTerm.(μv, ΔRR₀, sameRR₀) for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) uvwSum = u + v + w γ = μSum - uvwSum - term2Coeff2 = termProd5(uvwSum, T) * 2Fγs[γ+1] * - α^(rstSum - opqSum - uvwSum) - @inbounds tmp += mapMapReduce((u, v, w), core2s) * - term2Coeff1 * term2Coeff2 + @inbounds term2Coeff = subterm3 * termProd5(uvwSum, T) * 2Fγs[γ+1] * + α^(rstSum - opqSum - uvwSum) + term2 += mapMapReduce((u, v, w), genTerm2s) * term2Coeff end - A += mapMapReduce(rst, core1s) * tmp * - (term1Coeff2 * term1Coeff4 / (term1Coeff1 * term1Coeff3)) + A += mapMapReduce(rst, genTerm1s) * term1Coeff * term2 end end @@ -302,8 +301,8 @@ function ∫eeInteractionCore1234(ΔRl::NTuple{3, T}, ΔRr::NTuple{3, T}, lmn₄Sum = sum(lmn₄) lmnLSum = lmn₁Sum + lmn₂Sum lmnRSum = lmn₃Sum + lmn₄Sum - term1Coeff1 = termProd1(lmn₁, lmn₂) - term2Coeff1 = termProd1(lmn₄, lmn₃) + subterm1L = termProd1(lmn₁, lmn₂) + subterm1R = termProd1(lmn₄, lmn₃) for o₁ in 0:(i₁-2l₁), p₁ in 0:(j₁-2m₁), q₁ in 0:(k₁-2n₁), o₂ in 0:(i₂-2l₂), p₂ in 0:(j₂-2m₂), q₂ in 0:(k₂-2n₂), @@ -320,21 +319,19 @@ function ∫eeInteractionCore1234(ΔRl::NTuple{3, T}, ΔRr::NTuple{3, T}, opq₂Sum = sum(opq₂) opq₃Sum = sum(opq₃) opq₄Sum = sum(opq₄) - term1Coeff2 = termProd2(opqL) - term2Coeff2 = termProd2(opqR) - term1Coeff3 = termProd3(ijk₁, lmn₁, opq₁, ijk₂, lmn₂, opq₂) - term2Coeff3 = termProd3(ijk₄, lmn₄, opq₄, ijk₃, lmn₃, opq₃) + subterm2L = termProd3(ijk₁, lmn₁, opq₁, ijk₂, lmn₂, opq₂, T) / subterm1L + subterm2R = termProd3(ijk₄, lmn₄, opq₄, ijk₃, lmn₃, opq₃, T) / subterm1R μˣ, μʸ, μᶻ = μv = begin @. IJK - muladd(2, lmn₁+lmn₂+lmn₃+lmn₄, opqL+opqR) end - term3Coeff1 = termProd2(μv) + subterm3 = termProd2(μv) μsum = sum(μv) Fγs = F₀toFγ(μsum, β) - core1s = genIntTerm.(opqL, ΔRl, sameRl) - core2s = genIntTerm.(opqR, ΔRr, sameRr) - core3s = genIntTerm.(μv, ΔRc, sameRc) + genTerm1s = genIntTerm.(opqL, ΔRl, sameRl) + genTerm2s = genIntTerm.(opqR, ΔRr, sameRr) + genTerm3s = genIntTerm.(μv, ΔRc, sameRc) for r₁ in 0:((o₁+o₂)÷2), s₁ in 0:((p₁+p₂)÷2), t₁ in 0:((q₁+q₂)÷2), r₂ in 0:((o₃+o₄)÷2), s₂ in 0:((p₃+p₄)÷2), t₂ in 0:((q₃+q₄)÷2) @@ -343,21 +340,23 @@ function ∫eeInteractionCore1234(ΔRl::NTuple{3, T}, ΔRr::NTuple{3, T}, rst₂ = (r₂, s₂, t₂) rst₁Sum = sum(rst₁) rst₂Sum = sum(rst₂) - term1Coeff4 = termProd4(lmn₁Sum, opq₁Sum, lmn₂Sum, opq₂Sum, rst₁Sum, (α₁, α₂)) * (αL)^muladd(2, lmnLSum, rst₁Sum) - term2Coeff4 = termProd4(lmn₄Sum, opq₄Sum, lmn₃Sum, opq₃Sum, rst₂Sum, (α₄, α₃)) * (αR)^muladd(2, lmnRSum, rst₂Sum) - tmp = T(0.0) + term1Coeff = termProd4(lmn₁Sum, opq₁Sum, lmn₂Sum, opq₂Sum, rst₁Sum, + (α₁, α₂)) * (αL)^muladd(2, lmnLSum, rst₁Sum) * + subterm2L + term2Coeff = termProd4(lmn₄Sum, opq₄Sum, lmn₃Sum, opq₃Sum, rst₂Sum, + (α₄, α₃)) * (αR)^muladd(2, lmnRSum, rst₂Sum) * + subterm2R + term3 = T(0.0) for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) uvwSum = u + v + w γ = μsum - uvwSum - term3Coeff2 = termProd5(uvwSum, T) * 2Fγs[γ+1] * η^γ - @inbounds tmp += mapMapReduce((u, v, w), core3s) * - term3Coeff1 * term3Coeff2 + @inbounds term3Coeff = subterm3 * termProd5(uvwSum, T) * 2Fγs[γ+1] * η^γ + term3 += mapMapReduce((u, v, w), genTerm3s) * term3Coeff end - A += mapMapReduce(rst₁, core1s) * mapMapReduce(rst₂, core2s) * - (term1Coeff2 * term1Coeff4 / (term1Coeff1 * term1Coeff3)) * - (term2Coeff2 * term2Coeff4 / (term2Coeff1 * term2Coeff3)) * - tmp + A += mapMapReduce(rst₁, genTerm1s) * term1Coeff * + mapMapReduce(rst₂, genTerm2s) * term2Coeff * + term3 end end From adcb2b572fbc6f91bb73b68e648b1656824447f9 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Sun, 17 Mar 2024 19:04:43 -0400 Subject: [PATCH 18/22] Reduced useless compuation when a term is zero. --- src/Integrals/Core.jl | 90 +++++++++++++++++++++++++++---------------- 1 file changed, 57 insertions(+), 33 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index c7676fa7..7bceec47 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -226,29 +226,41 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T subterm3 = termProd2(μv) μSum = sum(μv) Fγs = @inbounds Fγss[μSum+1] + genTerm1s = genIntTerm.(opq₁.+opq₂, ΔR₁R₂, sameR₁R₂) + genTerm2s = genIntTerm.(μv, ΔRR₀, sameRR₀ ) for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) rst = (r, s, t) - rstSum = sum(rst) - term1Coeff = termProd4(lmn₁Sum, opq₁Sum, lmn₂Sum, opq₂Sum, rstSum, - (α₁, α₂)) * subterm2 - term2 = T(0.0) - genTerm2s = genIntTerm.(μv, ΔRR₀, sameRR₀) - - for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) - uvwSum = u + v + w - γ = μSum - uvwSum - @inbounds term2Coeff = subterm3 * termProd5(uvwSum, T) * 2Fγs[γ+1] * - α^(rstSum - opqSum - uvwSum) - term2 += mapMapReduce((u, v, w), genTerm2s) * term2Coeff + term1base = mapMapReduce(rst, genTerm1s) + + if !iszero(term1base) + + rstSum = sum(rst) + term1Coeff = termProd4(lmn₁Sum, opq₁Sum, lmn₂Sum, opq₂Sum, rstSum, + (α₁, α₂)) * subterm2 + term2 = T(0.0) + + for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) + + uvwSum = u + v + w + γ = μSum - uvwSum + term2base = mapMapReduce((u, v, w), genTerm2s) + + if !iszero(term2base) + @inbounds term2Coeff = subterm3 * termProd5(uvwSum, T) * + α^(rstSum - opqSum - uvwSum) * 2Fγs[γ+1] + term2 += term2base * term2Coeff + end + end + + A += term1base * term1Coeff * term2 end - A += mapMapReduce(rst, genTerm1s) * term1Coeff * term2 end end - end + A end @@ -338,29 +350,41 @@ function ∫eeInteractionCore1234(ΔRl::NTuple{3, T}, ΔRr::NTuple{3, T}, rst₁ = (r₁, s₁, t₁) rst₂ = (r₂, s₂, t₂) - rst₁Sum = sum(rst₁) - rst₂Sum = sum(rst₂) - term1Coeff = termProd4(lmn₁Sum, opq₁Sum, lmn₂Sum, opq₂Sum, rst₁Sum, - (α₁, α₂)) * (αL)^muladd(2, lmnLSum, rst₁Sum) * - subterm2L - term2Coeff = termProd4(lmn₄Sum, opq₄Sum, lmn₃Sum, opq₃Sum, rst₂Sum, - (α₄, α₃)) * (αR)^muladd(2, lmnRSum, rst₂Sum) * - subterm2R - term3 = T(0.0) - - for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) - uvwSum = u + v + w - γ = μsum - uvwSum - @inbounds term3Coeff = subterm3 * termProd5(uvwSum, T) * 2Fγs[γ+1] * η^γ - term3 += mapMapReduce((u, v, w), genTerm3s) * term3Coeff + term1base = mapMapReduce(rst₁, genTerm1s) + term2base = mapMapReduce(rst₂, genTerm2s) + term12base = term1base * term2base + + if !iszero(term12base) + + rst₁Sum = sum(rst₁) + rst₂Sum = sum(rst₂) + term1Coeff = termProd4(lmn₁Sum, opq₁Sum, lmn₂Sum, opq₂Sum, rst₁Sum, + (α₁, α₂)) * (αL)^muladd(2, lmnLSum, rst₁Sum) * + subterm2L + term2Coeff = termProd4(lmn₄Sum, opq₄Sum, lmn₃Sum, opq₃Sum, rst₂Sum, + (α₄, α₃)) * (αR)^muladd(2, lmnRSum, rst₂Sum) * + subterm2R + term3 = T(0.0) + + for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) + + uvwSum = u + v + w + γ = μsum - uvwSum + term3base = mapMapReduce((u, v, w), genTerm3s) + + if !iszero(term3base) + @inbounds term3Coeff = subterm3 * termProd5(uvwSum, T) * + η^γ * 2Fγs[γ+1] + term3 += term3base * term3Coeff + end + end + + A += term12base * term1Coeff * term2Coeff * term3 end - A += mapMapReduce(rst₁, genTerm1s) * term1Coeff * - mapMapReduce(rst₂, genTerm2s) * term2Coeff * - term3 end end - end + A end From 1dd74a42c276a004b62ee9739525f2fadad3fb52 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 18 Mar 2024 03:06:32 -0400 Subject: [PATCH 19/22] Attempt for better performance based on `getIntRange`. --- src/Integrals/Core.jl | 135 ++++++++++++++++++++++++++---------------- 1 file changed, 84 insertions(+), 51 deletions(-) diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 7bceec47..9a682e01 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -67,39 +67,52 @@ function F₀toFγ(γ::Int, u::T, resHolder::Vector{T}=Array{T}(undef, γ+1)) wh resHolder end +## Not very efficient. +# function getIntRange(doubleUpper::T1, zeroΔx::Bool) where {T1<:Integer} +# if zeroΔx +# if iseven(doubleUpper) +# let Ω=doubleUpper÷2 +# Ω:Ω +# end +# else +# T1(1):T1(0) +# end +# else +# T1(0):(doubleUpper÷2) +# end +# end + +# function genIntTerm(doubleUpper::T1, Δx::T2, zeroΔx::Bool) where {T1<:Integer, T2<:Real} +# f = if zeroΔx +# function (arg::T1) +# (inv∘T2∘factorial)(arg) +# end +# else +# let doubleUpper=doubleUpper, Δx=Δx +# function (arg::T1) +# pwr = muladd(-2, arg, doubleUpper) +# Δx^pwr / factorial(pwr) / factorial(arg) +# end +# end +# end +# f +# end -function genIntOverlapTerm(getRangeFunc::F, Δx::T, i₁::Int, α₁::T, - i₂::Int, α₂::T) where {F, T} - f = let getRange=getRangeFunc, Δx=Δx, i₁=i₁, α₁=α₁, i₂=i₂, α₂=α₂ - function (l₁::Int, l₂::Int) - res = T(0.0) - Ω = muladd(-2, l₁ + l₂, i₁ + i₂) - for o in getRange(Ω) - res += Δx^(Ω-2o) * - ( α₁^(i₂ - l₁ - 2l₂ - o) / (factorial(l₂) * factorial(i₁-2l₁)) ) * - ( α₂^(i₁ - l₂ - 2l₁ - o) / (factorial(l₁) * factorial(i₂-2l₂)) ) * - ( n1Power(o) * factorial(Ω) * (α₁ + α₂)^muladd(2, (l₁ + l₂), o) / - ((1 << 2(l₁ + l₂ + o)) * factorial(o ) * factorial(Ω -2o )) ) +function genIntTerm(arg1::T1, Δx::T2, zeroΔx::Bool) where {T1<:Integer, T2<:Real} + fRes = let arg1=arg1, Δx=Δx, zeroΔx=zeroΔx + function (arg2::T1) + pwr = muladd(-2, arg2, arg1) + coeff = if zeroΔx + iszero(pwr) ? T2(1.0) : (return T2(0.0)) + else + Δx^pwr / factorial(pwr) end - res + coeff / factorial(arg2) end end - f + fRes end -function genIntOverlapCore(Δx::T, i₁::Int, α₁::T, - i₂::Int, α₂::T, - zeroΔx::Bool=false) where {T} - getRange = ifelse( zeroΔx, - Ω::Int -> ifelse(iseven(Ω), begin ΩHalf=Ω÷2; ΩHalf:ΩHalf end, 1:0), - Ω::Int -> 0:(Ω÷2) ) - f = genIntOverlapTerm(getRange, Δx, i₁, α₁, i₂, α₂) - res = T(0.0) - for l₁ in 0:(i₁÷2), l₂ in 0:(i₂÷2) - res += f(l₁, l₂) - end - res -end function ∫overlapCore(::Val{3}, ΔR::NTuple{3, T}, @@ -109,14 +122,44 @@ function ∫overlapCore(::Val{3}, coeff::T=T(1.0)) where {T} any(n -> n<0, (ijk₁..., ijk₂...)) && (return T(0.0)) α = α₁ + α₂ - res = T(1.0) - for (i₁, i₂, ΔRᵢ, zeroΔx) in zip(ijk₁, ijk₂, ΔR, sameCen) - int = genIntOverlapCore(ΔRᵢ, i₁, α₁, i₂, α₂, zeroΔx) - iszero(int) && (return T(0.0)) - res *= n1Power(i₁) * factorial(i₁) * factorial(i₂) * α^(-i₁-i₂) * int + int = T(1.0) + + for (i₁, i₂, Δx, zeroΔx) in zip(ijk₁, ijk₂, ΔR, sameCen) + + i = i₁ + i₂ + res = T(0.0) + + for l₁ in 0:(i₁÷2), l₂ in 0:(i₂÷2) + + l = l₁ + l₂ + Ω = muladd(-2, l, i) + subterm1 = (T∘factorial)(Ω) / ( factorial(l₁) * factorial(i₁-2l₁) * + factorial(l₂) * factorial(i₂-2l₂) ) + term = T(0.0) + # oRange = getIntRange(Ω, zeroΔx) + genTerm = genIntTerm(Ω, Δx, zeroΔx) + + for o in 0:(Ω÷2) + termBase = genTerm(o) + if !iszero(termBase) + subterm2 = n1Power(o) * (α₁ + α₂)^muladd(2, l, o) * + α₁^(i₂ - l - l₂ - o) * α₂^(i₁ - l - l₁ - o) / (1 << 2(l + o)) + term += termBase * subterm1 * subterm2 + end + end + + res += term + end + + if iszero(res) + return T(0.0) + else + int *= res + end end - res *= sqrt((π/α)^3) * exp(-α₁ / α * α₂ * sum(abs2, ΔR)) - res * coeff + + int * prod(@. n1Power(ijk₁) * factorial(ijk₁) * factorial(ijk₂) / α^(ijk₁ + ijk₂)) * + sqrt((π/α)^3) * exp(-α₁ / α * α₂ * sum(abs2, ΔR)) * coeff end ∫overlapCore(::Val{3}, @@ -149,21 +192,6 @@ function ∫elecKineticCore(::Val{3}, end -function genIntTerm(arg1::T1, Δx::T2, zeroΔx::Bool) where {T1<:Integer, T2<:Real} - fRes = let arg1=arg1, Δx=Δx, zeroΔx=zeroΔx - function (arg2::T1) - pwr = muladd(-2, arg2, arg1) - coeff = if zeroΔx - iszero(pwr) ? T2(1.0) : (return T2(0.0)) - else - Δx^pwr / factorial(pwr) - end - coeff / factorial(arg2) - end - end - fRes -end - function termProd1(lmn₁::NTuple{3, T1}, lmn₂::NTuple{3, T1}) where {T1<:Integer} prod( @. factorial(lmn₁) * factorial(lmn₂) ) end @@ -176,8 +204,8 @@ function termProd3(ijk₁::NTuple{3, T1}, lmn₁::NTuple{3, T1}, opq₁::NTuple{ ijk₂::NTuple{3, T1}, lmn₂::NTuple{3, T1}, opq₂::NTuple{3, T1}, ::Type{T2}) where {T1<:Integer, T2<:Real} (T2∘termProd2)(opq₁ .+ opq₂) / - prod( @. ( factorial(opq₁) * factorial(ijk₁-2lmn₁-opq₁) * - factorial(opq₂) * factorial(ijk₂-2lmn₂-opq₂) ) ) + prod( @. factorial(opq₁) * factorial(ijk₁-2lmn₁-opq₁) * + factorial(opq₂) * factorial(ijk₂-2lmn₂-opq₂) ) end function termProd4(lmn₁Sum::T1, opq₁Sum::T1, lmn₂Sum::T1, opq₂Sum::T1, rstSum::T1, (α₁, α₂)::NTuple{2, T2}) where {T1<:Integer, T2<:Real} @@ -222,14 +250,18 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T opqSum = opq₁Sum + opq₂Sum subterm2 = termProd3(ijk₁, lmn₁, opq₁, ijk₂, lmn₂, opq₂, T) / subterm1 + # μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq) μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq) subterm3 = termProd2(μv) μSum = sum(μv) Fγs = @inbounds Fγss[μSum+1] + # rstRange = getIntRange.(opq₁.+opq₂, sameR₁R₂) genTerm1s = genIntTerm.(opq₁.+opq₂, ΔR₁R₂, sameR₁R₂) + # uvwRange = getIntRange.(μv, sameRR₀ ) genTerm2s = genIntTerm.(μv, ΔRR₀, sameRR₀ ) + # for r in rstRange[1], s in rstRange[2], t in rstRange[3] for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) rst = (r, s, t) @@ -242,6 +274,7 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T (α₁, α₂)) * subterm2 term2 = T(0.0) + # for u in uvwRange[1], v in uvwRange[2], w in uvwRange[3] for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) uvwSum = u + v + w From 46cfd304f2202e284a7c77d8b01a0e8c0b1029fa Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Mon, 18 Mar 2024 16:06:24 -0400 Subject: [PATCH 20/22] Code reformatting. --- src/Basis.jl | 2 +- src/Integrals/Core.jl | 111 ++++++++++++++++-------------------------- src/Library.jl | 3 +- src/Tools.jl | 4 +- 4 files changed, 49 insertions(+), 71 deletions(-) diff --git a/src/Basis.jl b/src/Basis.jl index 8b78c6e3..9cb92b93 100644 --- a/src/Basis.jl +++ b/src/Basis.jl @@ -2070,7 +2070,7 @@ hasNormFactor(b::FloatingGTBasisFuncs) = b.normalizeGTO getNijk(::Type{T}, i::Integer, j::Integer, k::Integer) where {T} = -T(πvals[-0.75]) * 2^T(1.5*(i+j+k) + 0.75) * T(sqrt( factorial(i) * factorial(j) * +T(πPowers[:n0d75]) * 2^T(1.5*(i+j+k) + 0.75) * T(sqrt( factorial(i) * factorial(j) * factorial(k) / (factorial(2i) * factorial(2j) * factorial(2k)) )) getNα(i::Integer, j::Integer, k::Integer, α::T) where {T} = diff --git a/src/Integrals/Core.jl b/src/Integrals/Core.jl index 9a682e01..0e0d4a42 100644 --- a/src/Integrals/Core.jl +++ b/src/Integrals/Core.jl @@ -40,7 +40,7 @@ function F0(u::T) where {T} T(1.0) else ur = sqrt(u) - T(πvals[0.5]) * erf(ur) / (2ur) + T(πPowers[:p0d5]) * erf(ur) / (2ur) end end @@ -67,36 +67,6 @@ function F₀toFγ(γ::Int, u::T, resHolder::Vector{T}=Array{T}(undef, γ+1)) wh resHolder end -## Not very efficient. -# function getIntRange(doubleUpper::T1, zeroΔx::Bool) where {T1<:Integer} -# if zeroΔx -# if iseven(doubleUpper) -# let Ω=doubleUpper÷2 -# Ω:Ω -# end -# else -# T1(1):T1(0) -# end -# else -# T1(0):(doubleUpper÷2) -# end -# end - -# function genIntTerm(doubleUpper::T1, Δx::T2, zeroΔx::Bool) where {T1<:Integer, T2<:Real} -# f = if zeroΔx -# function (arg::T1) -# (inv∘T2∘factorial)(arg) -# end -# else -# let doubleUpper=doubleUpper, Δx=Δx -# function (arg::T1) -# pwr = muladd(-2, arg, doubleUpper) -# Δx^pwr / factorial(pwr) / factorial(arg) -# end -# end -# end -# f -# end function genIntTerm(arg1::T1, Δx::T2, zeroΔx::Bool) where {T1<:Integer, T2<:Real} fRes = let arg1=arg1, Δx=Δx, zeroΔx=zeroΔx @@ -114,12 +84,19 @@ function genIntTerm(arg1::T1, Δx::T2, zeroΔx::Bool) where {T1<:Integer, T2<:Re end +computeUnit1(::Val{S}, divider::T) where {S, T} = T(πPowers[S]) / divider + +computeUnit2(η::T, ΔR::NTuple{N, T}) where {T, N} = exp(-η * sum(abs2, ΔR)) + +computeUnit3(ijk₁::NTuple{N, T1}, ijk₂::NTuple{N, T1}, α::T2) where {N, T1, T2} = +prod(@. factorial(ijk₁) * factorial(ijk₂) / α^(ijk₁ + ijk₂)) + + function ∫overlapCore(::Val{3}, ΔR::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - sameCen::NTuple{3, Bool}=(false, false, false), - coeff::T=T(1.0)) where {T} + sameCen::NTuple{3, Bool}=(false, false, false)) where {T} any(n -> n<0, (ijk₁..., ijk₂...)) && (return T(0.0)) α = α₁ + α₂ int = T(1.0) @@ -136,14 +113,13 @@ function ∫overlapCore(::Val{3}, subterm1 = (T∘factorial)(Ω) / ( factorial(l₁) * factorial(i₁-2l₁) * factorial(l₂) * factorial(i₂-2l₂) ) term = T(0.0) - # oRange = getIntRange(Ω, zeroΔx) genTerm = genIntTerm(Ω, Δx, zeroΔx) for o in 0:(Ω÷2) termBase = genTerm(o) if !iszero(termBase) subterm2 = n1Power(o) * (α₁ + α₂)^muladd(2, l, o) * - α₁^(i₂ - l - l₂ - o) * α₂^(i₁ - l - l₁ - o) / (1 << 2(l + o)) + α₁^(i₂ - l - l₂ - o) * α₂^(i₁ - l - l₁ - o) / (1 << 2(l + o)) term += termBase * subterm1 * subterm2 end end @@ -158,17 +134,17 @@ function ∫overlapCore(::Val{3}, end end - int * prod(@. n1Power(ijk₁) * factorial(ijk₁) * factorial(ijk₂) / α^(ijk₁ + ijk₂)) * - sqrt((π/α)^3) * exp(-α₁ / α * α₂ * sum(abs2, ΔR)) * coeff + η = α₁ / α * α₂ + int * computeUnit1(Val(:p1d5), α^T(1.5)) * computeUnit2(η, ΔR) * + computeUnit3(ijk₁, ijk₂, α) * n1Power(ijk₁) end ∫overlapCore(::Val{3}, R₁::NTuple{3, T}, R₂::NTuple{3, T}, ijk₁::NTuple{3, Int}, α₁::T, ijk₂::NTuple{3, Int}, α₂::T, - sameCen::NTuple{3, Bool}=(false, false, false), - coeff::T=T(1.0)) where {T} = -∫overlapCore(Val(3), R₁.-R₂, ijk₁, α₁, ijk₂, α₂, sameCen, coeff) + sameCen::NTuple{3, Bool}=(false, false, false)) where {T} = +∫overlapCore(Val(3), R₁.-R₂, ijk₁, α₁, ijk₂, α₂, sameCen) function ∫elecKineticCore(::Val{3}, @@ -178,17 +154,17 @@ function ∫elecKineticCore(::Val{3}, sameCen::NTuple{3, Bool}=(false, false, false)) where {T} ΔR = R₁ .- R₂ shifts = ((1,0,0), (0,1,0), (0,0,1)) - mapreduce(+, ijk₁, ijk₂, shifts) do i₁, i₂, Δ𝑙 + map(ijk₁, ijk₂, shifts) do i₁, i₂, Δ𝑙 Δijk1 = map(-, ijk₁, Δ𝑙) Δijk2 = map(-, ijk₂, Δ𝑙) Δijk3 = map(+, ijk₁, Δ𝑙) Δijk4 = map(+, ijk₂, Δ𝑙) - int1 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk2, α₂, sameCen, T(i₁*i₂/2)) - int2 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk4, α₂, sameCen, 2α₁*α₂ ) - int3 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk2, α₂, sameCen, α₁*i₂ ) - int4 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk4, α₂, sameCen, α₂*i₁ ) - int1 + int2 - int3 - int4 - end + int1 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk2, α₂, sameCen) + int2 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk4, α₂, sameCen) + int3 = ∫overlapCore(Val(3), ΔR, Δijk3, α₁, Δijk2, α₂, sameCen) + int4 = ∫overlapCore(Val(3), ΔR, Δijk1, α₁, Δijk4, α₂, sameCen) + int1*i₁*i₂/2 + int2*2α₁*α₂ - int3*α₁*i₂ - int4*α₂*i₁ + end |> sum end @@ -250,18 +226,14 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T opqSum = opq₁Sum + opq₂Sum subterm2 = termProd3(ijk₁, lmn₁, opq₁, ijk₂, lmn₂, opq₂, T) / subterm1 - # μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq) μˣ, μʸ, μᶻ = μv = @. ijk₁ + ijk₂ - muladd(2, lmn₁+lmn₂, opq) subterm3 = termProd2(μv) μSum = sum(μv) Fγs = @inbounds Fγss[μSum+1] - # rstRange = getIntRange.(opq₁.+opq₂, sameR₁R₂) genTerm1s = genIntTerm.(opq₁.+opq₂, ΔR₁R₂, sameR₁R₂) - # uvwRange = getIntRange.(μv, sameRR₀ ) genTerm2s = genIntTerm.(μv, ΔRR₀, sameRR₀ ) - # for r in rstRange[1], s in rstRange[2], t in rstRange[3] for r in 0:((o₁+o₂)÷2), s in 0:((p₁+p₂)÷2), t in 0:((q₁+q₂)÷2) rst = (r, s, t) @@ -274,7 +246,6 @@ function genIntNucAttCore(ΔRR₀::NTuple{3, T}, ΔR₁R₂::NTuple{3, T}, β::T (α₁, α₂)) * subterm2 term2 = T(0.0) - # for u in uvwRange[1], v in uvwRange[2], w in uvwRange[3] for u in 0:(μˣ÷2), v in 0:(μʸ÷2), w in 0:(μᶻ÷2) uvwSum = u + v + w @@ -305,14 +276,17 @@ function ∫nucAttractionCore(::Val{3}, sameR₁R₂::NTuple{3, Bool}=(false, false, false)) where {T} α = α₁ + α₂ R = @. (α₁*R₁ + α₂*R₂) / α - ΔRR₀ = R .- R₀ + ΔRR₀ = R .- R₀ ΔR₁R₂ = R₁ .- R₂ β = α * sum(abs2, ΔRR₀) sameRR₀ = iszero.(ΔRR₀) - genIntNucAttCore(ΔRR₀, ΔR₁R₂, β, ijk₁, α₁, ijk₂, α₂, (sameRR₀, sameR₁R₂)) * - (π / α) * exp(-α₁ / α * α₂ * sum(abs2, ΔR₁R₂)) * - ( -Z₀ * (-1)^sum(ijk₁ .+ ijk₂) * - mapMapReduce(ijk₁, factorial) * mapMapReduce(ijk₂, factorial) ) + res = genIntNucAttCore(ΔRR₀, ΔR₁R₂, β, ijk₁, α₁, ijk₂, α₂, (sameRR₀, sameR₁R₂)) + if !iszero(res) + η = α₁ / α * α₂ + res *= computeUnit1(Val(:p1d0), α) * computeUnit2(η, ΔR₁R₂) * + prod(@. factorial(ijk₁) * factorial(ijk₂)) * (-Z₀ * n1Power(ijk₁ .+ ijk₂)) + end + res end @@ -429,24 +403,25 @@ function ∫eeInteractionCore(::Val{3}, R₄::NTuple{3, T}, ijk₄::NTuple{3, Int}, α₄::T, (sameR₁R₂, sameR₃R₄)::NTuple{2, NTuple{3, Bool}}= (Tuple∘fill)((false, false, false), 2)) where {T} - ΔRl = R₁ .- R₂ - ΔRr = R₃ .- R₄ αl = α₁ + α₂ αr = α₃ + α₄ - ηl = α₁ / αl * α₂ - ηr = α₃ / αr * α₄ + ΔRl = R₁ .- R₂ + ΔRr = R₃ .- R₄ ΔRc = @. (α₁*R₁ + α₂*R₂) / αl - (α₃*R₃ + α₄*R₄) / αr zeroΔRc = iszero.(ΔRc) η = αl / (α₁ + α₂ + α₃ + α₄) * αr β = η * sum(abs2, ΔRc) - ∫eeInteractionCore1234(ΔRl, ΔRr, ΔRc, β, η, ijk₁, α₁, ijk₂, α₂, ijk₃, α₃, ijk₄, α₄, - (sameR₁R₂, sameR₃R₄, zeroΔRc)) * - T(πvals[2.5]) / (αl * αr * sqrt(αl + αr)) * - exp(-ηl * sum(abs2, ΔRl)) * exp(-ηr * sum(abs2, ΔRr)) * - mapreduce(*, ijk₁, ijk₂, ijk₃, ijk₄) do l₁, l₂, l₃, l₄ - (-1)^(l₁+l₂) * factorial(l₁) * factorial(l₂) * factorial(l₃) * factorial(l₄) / - (αl^(l₁+l₂) * αr^(l₃+l₄)) + res = ∫eeInteractionCore1234(ΔRl, ΔRr, ΔRc, β, η, + ijk₁, α₁, ijk₂, α₂, ijk₃, α₃, ijk₄, α₄, + (sameR₁R₂, sameR₃R₄, zeroΔRc)) + if !iszero(res) + ηl = α₁ / αl * α₂ + ηr = α₃ / αr * α₄ + res *= computeUnit1(Val(:p2d5), αl*αr*sqrt(αl + αr)) * + computeUnit2(ηl, ΔRl) * computeUnit3(ijk₁, ijk₂, αl) * + computeUnit2(ηr, ΔRr) * computeUnit3(ijk₃, ijk₄, αr) * n1Power(ijk₁ .+ ijk₂) end + res end diff --git a/src/Library.jl b/src/Library.jl index c9e2347f..991e2afb 100644 --- a/src/Library.jl +++ b/src/Library.jl @@ -247,7 +247,8 @@ function checkBSList(;printInfo::Bool=false) end -const πvals = Dict([-0.75, 0.5, 2.5] .=> big(π).^[-0.75, 0.5, 2.5]) +const πPowers = Dict( [:n0d75, :p0d5, :p1d0, :p1d5, :p2d5] .=> big(π).^ + [ -0.75, 0.5, 1.0, 1.5, 2.5] ) const DefaultDigits = 10 diff --git a/src/Tools.jl b/src/Tools.jl index 7df942ed..857f3d68 100644 --- a/src/Tools.jl +++ b/src/Tools.jl @@ -1030,4 +1030,6 @@ function betterSum(x) # Improved Kahan–Babuška algorithm fro array summation end n1Power(x::Int) = ifelse(isodd(x), -1, 1) -n1Power(x::Int, ::Type{T}) where {T} = ifelse(isodd(x), T(-1), T(1)) \ No newline at end of file +n1Power(x::Int, ::Type{T}) where {T} = ifelse(isodd(x), T(-1), T(1)) +n1Power(x::NTuple{N, Int}) where {N} = (n1Power∘sum)(x) +n1Power(x::NTuple{N, Int}, ::Type{T}) where {N, T} = n1Power(sum(x), T) \ No newline at end of file From 331632bfbfdfda5a91f7f09311ad1cb073e89de1 Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Wed, 20 Mar 2024 15:55:08 -0400 Subject: [PATCH 21/22] Merged test dependency into the main Project.toml. --- Project.toml | 21 +++++++++++++++++++++ test/Project.toml | 13 ------------- 2 files changed, 21 insertions(+), 13 deletions(-) delete mode 100644 test/Project.toml diff --git a/Project.toml b/Project.toml index 02fe0e25..a9cebaa7 100644 --- a/Project.toml +++ b/Project.toml @@ -20,6 +20,7 @@ TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" [compat] Combinatorics = "1" +Documenter = "1" DoubleFloats = "^1.2" FastGaussQuadrature = "1" ForwardDiff = "^0.10.25" @@ -28,9 +29,29 @@ LRUCache = "^1.6" LazyArrays = "^1.8" LineSearches = "^7.1.1" LinearAlgebra = "1" +Optim = "1" Printf = "1" +QuadGK = "2" +Random = "1" SPGBox = "0.7" SpecialFunctions = "2" Statistics = "1" +Suppressor = "0.2" TensorOperations = "^4.0.5" +Test = "1" julia = "^1.6" +libcint_jll = "5" + +[extras] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +libcint_jll = "574b78ca-bebd-517c-801d-4735c93a9686" + +[targets] +test = ["Documenter", "ForwardDiff", "LinearAlgebra", "Optim", "QuadGK", "Random", "Suppressor", "Test", "libcint_jll"] \ No newline at end of file diff --git a/test/Project.toml b/test/Project.toml deleted file mode 100644 index 71af9f5c..00000000 --- a/test/Project.toml +++ /dev/null @@ -1,13 +0,0 @@ -[deps] -Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Optim = "429524aa-4258-5aef-a3af-852621145aeb" -QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" -Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" -Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -libcint_jll = "574b78ca-bebd-517c-801d-4735c93a9686" - -[compat] -libcint_jll = "5" \ No newline at end of file From 0c872d88a12afb509d083e6eb2b2565b632f12fe Mon Sep 17 00:00:00 2001 From: frankwswang <41162655+frankwswang@users.noreply.github.com> Date: Wed, 20 Mar 2024 15:55:54 -0400 Subject: [PATCH 22/22] Updated the package version number. --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a9cebaa7..6da3ce45 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Quiqbox" uuid = "7cb8c394-fae1-4ab9-92f2-30189d7746cd" -version = "0.5.8" +version = "0.5.9" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"