Merge pull request #18 from ACEsuit/v0.1-dev

Cleanup towards v0.1

Author: cortner

.github/workflows/CI.yml

@@ -18,8 +18,8 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.9'
           - '1.10'
+          - '1.11'
         os:
           - ubuntu-latest
         arch:

.gitignore

@@ -1,2 +1,3 @@
 /Manifest.toml
 /docs/build/
+.vscode

Project.toml

@@ -1,7 +1,6 @@
 name = "RepLieGroups"
 uuid = "f07d36f2-91c4-427a-b67b-965fe5ebe1d2"
-authors = ["Christoph Ortner <christophortner@gmail.com> and contributors"]
-version = "0.0.4"
+version = "0.1.0-dev"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -13,19 +12,20 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 Combinatorics = "1.0"
-LinearAlgebra = "1.9"
+LinearAlgebra = "1.10"
 PartialWaveFunctions = "0.2.0"
-SparseArrays = "1.9"
+SparseArrays = "1.10"
 StaticArrays = "1.5"
-julia = "1.9, 1.10"
+julia = "1.10, 1.11"
 
 [extras]
 BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
-Polynomials4ML = "03c4bcba-a943-47e9-bfa1-b1661fc2974f"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 SpheriCart = "5caf2b29-02d9-47a3-9434-5931c85ba645"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 WignerD = "87c4ff3e-34df-11e9-37a7-516cea4e0402"
 
 [targets]
-test = ["Test", "Polynomials4ML", "WignerD", "Rotations", "BlockDiagonals", "SpheriCart"]
+test = ["Test", "WignerD", "Rotations", "BlockDiagonals", 
+        "SpheriCart", "Random"]

benchmarks/run_benchmark.jl

@@ -0,0 +1,96 @@
+
+using StaticArrays, LinearAlgebra, RepLieGroups, WignerD, Combinatorics, 
+      Rotations
+using WignerD: wignerD
+using Test
+
+isdefined(Main, :___UTILS_FOR_TESTS___) || include("../utils/utils_for_tests.jl")
+
+##
+
+# Test the new RPE basis up to L = 4
+@info("Testing the new cSH-based RPE basis")
+for L = 0:4
+   @info("Testing L = $L")
+
+   # generate a nnll_list for testing
+   lmax = 4
+   nmax = 4
+   nnll_list = [] 
+
+   for ORD = 2:6
+      for ll in with_replacement_combinations(1:lmax, ORD) 
+         # 0 or 1 above ?
+         if !iseven(sum(ll)+L); continue; end  # This is to ensure the reflection symmetry
+         if sum(ll) > 2 * lmax; continue; end 
+         for Inn in CartesianIndices( ntuple(_->1:nmax, ORD) )
+            nn = [ Inn.I[α] for α = 1:ORD ]
+            if sum(nn) > sum(1:nmax); continue; end
+            nnll = [ (ll[α], nn[α]) for α = 1:ORD ]
+            if !issorted(nnll); continue; end
+            push!(nnll_list, (SVector(nn...), SVector(ll...)))
+         end
+      end
+   end
+
+   long_nnll_list = nnll_list 
+   short_nnll_list = nnll_list[1:10:end]
+   ultra_short_nnll_list = nnll_list[1:100:end]
+
+   verbose = true 
+
+   @info("Using short nnll list for testing")
+   nnll_list = short_nnll_list
+
+   for (itest, (nn, ll)) in enumerate(nnll_list)
+      N = length(ll)
+      @assert length(ll) == length(nn)
+      
+      local Rs, θ, Q, QRs, D, D_r
+
+      Rs = rand_config(length(ll))
+      θ = rand(3) * 2pi
+      Q = RotZYZ(θ...)
+      D = transpose(wignerD(L, θ...)) 
+      D_r = Ctran(L) * D * Ctran(L)'
+      QRs = [Q*Rs[i] for i in 1:length(Rs)]
+      
+      t_rpe = @elapsed coeffs_ind, MM = O3.rpe_basis_new(nn,ll,L)
+      t_rpe_r = @elapsed coeffs_ind_r, MM_r = O3.rpe_basis_new(nn,ll,L; flag = :SpheriCart)
+      rk = rank(O3.gram(coeffs_ind),rtol = 1e-12)
+      rk_r = rank(O3.gram(coeffs_ind_r),rtol = 1e-12)
+      @test rk == rk_r
+
+      fRs1 = eval_basis(Rs; coeffs = coeffs_ind, MM = MM, ll = ll, nn = nn)
+      fRs1Q = eval_basis(QRs; coeffs = coeffs_ind, MM = MM, ll = ll, nn = nn)
+      # @info("Testing the equivariance of the old RPE basis")
+      L == 0 ? (@test norm(fRs1 - fRs1Q) < 1e-14) : (@test norm(fRs1 - Ref(D) .* fRs1Q) < 1e-14)
+
+      fRs1_r = eval_basis(Rs; coeffs = coeffs_ind_r, MM = MM_r, ll = ll, nn = nn, Real = true)
+      fRs1Q_r = eval_basis(QRs; coeffs = coeffs_ind_r, MM = MM_r, ll = ll, nn = nn, Real = true)
+      # @info("Testing the equivariance of the old RPE basis")
+      L == 0 ? (@test norm(fRs1_r - fRs1Q_r) < 1e-14) : (@test norm(fRs1_r - Ref(D_r) .* fRs1Q_r) < 1e-14)
+
+      ntest = 1000
+
+      RR = make_batch(ntest, length(ll))
+
+      X = rand_batch(; coeffs=coeffs_ind, MM=MM, ll=ll, nn=nn, batch = RR)
+      @test rank(O3.gram(X); rtol=1e-12) == size(X,1) == rk
+
+      X_r = rand_batch(; coeffs=coeffs_ind_r, MM=MM_r, ll=ll, nn=nn, batch = RR)
+      @test rank(O3.gram(X_r); rtol=1e-12) == size(X_r,1) == rk_r
+
+      Xsym = sym_rand_batch(; coeffs=coeffs_ind, MM=MM, ll=ll, nn=nn, batch = RR)
+      @test rank(O3.gram(Xsym); rtol=1e-12) == size(Xsym,1) == rk
+
+      Xsym_r = sym_rand_batch(; coeffs=coeffs_ind_r, MM=MM_r, ll=ll, nn=nn, batch = RR)
+      @test rank(O3.gram(Xsym_r); rtol=1e-12) == size(Xsym_r,1) == rk_r
+
+      if verbose 
+         @info("Test $itest: t_rpe = $t_rpe, t_rpe_r = $t_rpe_r")
+      else 
+         print(".")
+      end
+   end
+end

src/obsolete/O3.jl src/Archived/O3.jlsrc/obsolete/O3.jl renamed to src/Archived/O3.jl

@@ -257,7 +257,7 @@ function Ctran(i::Int64,j::Int64;convention = :SpheriCart)
 	end
 end
 
-Ctran(l::Int64;convention = :SpheriCart) = sparse(Matrix{ComplexF64}([ Ctran(m,μ;convention=convention) for m = -l:l, μ = -l:l ])) |> dropzeros
+Ctran(l::Int64; convention = :SpheriCart) = sparse(Matrix{ComplexF64}([ Ctran(m,μ;convention=convention) for m = -l:l, μ = -l:l ])) |> dropzeros
 
 ## NOTE: Ctran(L) is the transformation matrix from rSH to cSH. More specifically, 
 #        if we write Polynomials4ML rSH as R_{lm} and cSH as Y_{lm} and their corresponding 

src/obsolete/obsolete.jl src/Archived/obsolete.jlsrc/obsolete/obsolete.jl renamed to src/Archived/obsolete.jl

@@ -1,18 +1,54 @@
-module O3_new
-
-# Alternative to the computation of rotation equivariant coupling coefficients
+module O3
 
 import RepLieGroups
+
 using PartialWaveFunctions
 using Combinatorics
 using LinearAlgebra
 using StaticArrays
+using SparseArrays
+
+export coupling_coeffs
+
+
+# ------------------------------------------------------- 
+
+## NOTE: Ctran(L) is the transformation matrix from rSH to cSH. More specifically, 
+#        if we write Polynomials4ML rSH as R_{lm} and cSH as Y_{lm} and their corresponding 
+#        vectors of order L as R_L and Y_L, respectively. Then R_L = Ctran(L) * Y_L.
+#        This suggests that the "D-matrix" for the Polynomials4ML rSH is Ctran(l) * D(l) * Ctran(L)', 
+#        where D, the D-matrix for cSH. This inspires the following new CG recursion.
+
+# transformation matrix from RSH to CSH for different conventions
+function Ctran(i::Int64,j::Int64;convention = :SpheriCart)
+	if convention == :cSH
+		return i == j
+	end
+	
+	order_dict = Dict(:SpheriCart => [1,2,3,4], :CondonShortley => [4,3,2,1], :FHIaims => [4,2,3,1])
+	val_list = [(-1)^(i), im, (-1)^(i+1)*im, 1] ./ sqrt(2)
+	if abs(i) != abs(j)
+		return 0 
+	elseif i == j == 0
+		return 1
+	elseif i > 0 && j > 0
+		return val_list[order_dict[convention][1]]
+	elseif i < 0 && j < 0
+		return val_list[order_dict[convention][2]]
+	elseif i < 0 && j > 0
+		return val_list[order_dict[convention][3]]
+	elseif i > 0 && j < 0
+		return val_list[order_dict[convention][4]]
+	end
+end
+
+Ctran(l::Int64; convention = :SpheriCart) = sparse(Matrix{ComplexF64}([ Ctran(m,μ;convention=convention) for m = -l:l, μ = -l:l ])) |> dropzeros
 
-using RepLieGroups.O3: Ctran
 
-export re_rpe, rpe_basis_new
+# -----------------------------------------------------
 
-# The generalized Clebsch Gordan Coefficients; variables of this function are fully inherited from the first ACE paper
+# The generalized Clebsch Gordan Coefficients; variables of this function are 
+# fully inherited from the first ACE paper. 
 function GCG(l::SVector{N,Int64},m::SVector{N,Int64},L::SVector{N,Int64},M_N::Int64;flag=:cSH) where N
     # @assert -L[N] ≤ M_N ≤ L[N] 
     if m_filter(m, M_N;flag=flag) == false || L[1] < abs(m[1])
@@ -40,7 +76,7 @@ function GCG(l::SVector{N,Int64},m::SVector{N,Int64},L::SVector{N,Int64},M_N::In
                 mm = SA[mm...]
                 @assert sum(mm) == M
                 C_loc = GCG(l,mm,L,M;flag=:cSH)
-                coeff = Ctran(M_N,M;convention=flag)' * prod( RepLieGroups.O3.Ctran(m[i],mm[i];convention=flag) for i in 1:N )
+                coeff = Ctran(M_N,M;convention=flag)' * prod( Ctran(m[i],mm[i];convention=flag) for i in 1:N )
                 C_loc *= coeff
                 C += C_loc
             end
@@ -165,25 +201,113 @@ end
 # Function that generates the set of ordered m's given `n` and `l` with the abosolute sum of m's being smaller than L.
 m_generate(n,l,L;flag=:cSH) = union([m_generate(n,l,L,k;flag)[1] for k in -L:L]...), sum(length.(union([m_generate(n,l,L,k;flag)[1] for k in -L:L]...))) # orginal version: sum(m_generate(n,l,L,k;flag)[2] for k in -L:L), but this cannot be true anymore b.c. the m_classes can intersect
 
+function gram(X::Matrix{SVector{N,T}}) where {N,T}
+    G = zeros(T, size(X,1), size(X,1))
+    for i = 1:size(X,1)
+       for j = i:size(X,1)
+          G[i,j] = sum(dot(X[i,t], X[j,t]) for t = 1:size(X,2))
+          i == j ? nothing : (G[j,i]=G[i,j]')
+       end
+    end
+    return G
+ end
+
+gram(X::Matrix{<:Number}) = X * X'
+
+function lexi_ord(nn::SVector{N, Int64}, ll::SVector{N, Int64}) where N
+    pairs = collect(zip(ll, nn))         # create (lᵢ, nᵢ) pairs
+    sort!(pairs)                         # sort lexicographically: first by lᵢ, then by nᵢ
+    return SVector{N}(last.(pairs)), SVector{N}(first.(pairs))
+end
+
+"""
+    O3.coupling_coeffs(L, ll, nn; PI, basis)
+    O3.coupling_coeffs(L, ll; PI, basis)
+
+Compute coupling coefficients for the spherical harmonics basis, where 
+- `L` must be an `Integer`;
+- `ll, nn` must be vectors or tuples of `Integer` of the same length.
+- `PI`: whether or not the coupled basis is permutation-invariant (or the 
+corresponding tensor symmetric); default is `true` when `nn` is provided 
+and `false` when `nn` is not provided.
+- `basis`: which basis is being coupled, default is `complex`, alternative
+choice is `real`, which is compatible with the `SpheriCart.jl` convention.  
+"""
+function coupling_coeffs(L::Integer, ll, nn = nothing; 
+                         PI = !(isnothing(nn)), 
+                         basis = complex)
+
+    # convert L into the format required internally 
+    _L = Int(L) 
+
+    # convert ll into an SVector{N, Int}, as required internally 
+    N = length(ll) 
+    _ll = try 
+        _ll = SVector{N, Int}(ll...)
+    catch 
+        error("""coupling_coeffs(L::Integer, ll, ...) requires ll to be 
+               a vector or tuple of integers""")
+    end
 
-# Function that generates the coupling coefficient of the RE basis given `n` and `l`., the FMatrix for generating the RPE basis is also generated here.
-function re_rpe(n::SVector{N,Int64},l::SVector{N,Int64},L::Int64;flag=:cSH) where N
-    Lset = SetLl(l,L)
+    # convert nn into an SVector{N, Int}, as required internally 
+    if isnothing(nn) 
+        if PI 
+            _nn = SVector{N, Int}(ntuple(i -> 0, N)...)
+        else 
+            _nn = SVector{N, Int}((1:N)...)
+        end
+    elseif length(nn) != N 
+        error("""coupling_coeffs(L::Integer, ll, nn) requires ll and nn to be 
+               of the same length""")
+    else
+        _nn = try 
+            _nn = SVector{N, Int}(nn...)
+        catch 
+            error("""coupling_coeffs(L::Integer, ll, nn) requires nn to be 
+                   a vector or tuple of integers""")
+        end
+    end 
+
+    if basis == complex 
+        flag = :cSH 
+    elseif basis == real 
+        flag = :SpheriCart
+    elseif basis isa Symbol
+        flag = basis 
+    else 
+        error("unknown basis type: $basis")
+    end
+    
+    return _coupling_coeffs(_L, _ll, _nn; PI = PI, flag = flag)
+end
+    
+
+# Function that generates the coupling coefficient of the RE basis (PI = false) 
+# or RPE basis (PI = true) given `nn` and `ll`. 
+function _coupling_coeffs(L::Int64, ll::SVector{N, Int64}, nn::SVector{N, Int64}; 
+                          PI = true, flag = :cSH) where N
+
+    # TODO: when PI, (nn, ll) should be ordered 
+    if PI
+        nn, ll = lexi_ord(nn, ll)
+    end
+
+    Lset = SetLl(ll,L)
     r = length(Lset)
     T = L == 0 ? Float64 : SVector{2L+1,Float64}
     if r == 0 
         return zeros(T, 1, 1), zeros(T, 1, 1), [zeros(Int,N)], [zeros(Int,N)]
     else 
-        MMmat, size_m = m_generate(n,l,L;flag=flag) # classes of m's
+        MMmat, size_m = m_generate(nn,ll,L;flag=flag) # classes of m's
         FMatrix=zeros(T, r, length(MMmat)) # Matrix containing f(m,i)
         UMatrix=zeros(T, r, size_m) # Matrix containing the the coupling coefficients D
         MM = [] # all possible m's
         for i in 1:r
             c = 0
             for (j,m_class) in enumerate(MMmat)
-                for m in m_class
+                for mm in m_class
                     c += 1
-                    cg_coef = GCG(l,m,Lset[i];vectorize=(L!=0),flag=flag)
+                    cg_coef = GCG(ll,mm,Lset[i];vectorize=(L!=0),flag=flag)
                     FMatrix[i,j]+= cg_coef
                     UMatrix[i,c] = cg_coef
                 end
@@ -196,28 +320,16 @@ function re_rpe(n::SVector{N,Int64},l::SVector{N,Int64},L::Int64;flag=:cSH) wher
             end
         end      
     end
-    return UMatrix, FMatrix, MMmat, MM
-end
 
-function gram(X::Matrix{SVector{N,T}}) where {N,T}
-    G = zeros(T, size(X,1), size(X,1))
-    for i = 1:size(X,1)
-       for j = i:size(X,1)
-          G[i,j] = sum(dot(X[i,t], X[j,t]) for t = 1:size(X,2))
-          i == j ? nothing : (G[j,i]=G[i,j]')
-       end
+    if !PI
+        # return RE coupling coeffs if the permutation invariance is not needed
+        return UMatrix, MM
+    else
+        U, S, V = svd(gram(FMatrix))
+        rk = rank(Diagonal(S); rtol =  1e-12)
+        # return the RE-PI coupling coeffs
+        return Diagonal(S[1:rk]) * (U[:, 1:rk]' * UMatrix), MM
     end
-    return G
- end
-
- gram(X::Matrix{<:Number}) = X * X'
-
-function rpe_basis_new(nn::SVector{N, Int64}, ll::SVector{N, Int64}, L::Int64; flag = :cSH) where N
-    t_re = @elapsed UMatrix, FMatrix, MMmat, MM = re_rpe(nn, ll, L; flag = flag) # time of constructing the re_basis
-    # @show t_re # should be removed in the final version
-    U, S, V = svd(gram(FMatrix))
-    rk = rank(Diagonal(S); rtol =  1e-12)
-    return Diagonal(S[1:rk]) * (U[:, 1:rk]' * UMatrix), MM
 end
 
 end

test/side_tests.jl src/Archived/side_tests.jltest/side_tests.jl renamed to src/Archived/side_tests.jl

@@ -1,8 +1,8 @@
 module RepLieGroups
 
-include("yyvector.jl")
+export O3 
 
-include("obsolete/obsolete.jl")
+include("yyvector.jl")
 
 include("O3.jl")