Ms/equivarinat matrices #99

MatthiasSachs · 2022-01-05T15:26:25Z

All tests asserting the equivariance properties of EuclidianMatrix pass with tol = 1E-15.

The test println(@test(all(test_fio(basis; warntype = false)))) (see line 41 in test/test_EquivariantMatrix.jl) does NOT pass...

MatthiasSachs · 2022-01-05T15:46:10Z

test/test_EquivariantMatrix.jl not yet added to test collection. Thus, all tests are passing for now...

cortner · 2022-01-05T17:32:08Z

thank you - I'll take a look asap.

cortner · 2022-01-07T02:12:43Z

I'm looking into the read/write issue now.

cortner · 2022-01-07T02:13:24Z

Regarding the name, I think EquivariantMatrix is too general - it doesn't encode the specific symmetry. How about EuclideanMatrix? This would then also generalize nicely to EuclideanTensor if/when we get to that.

cortner · 2022-01-07T02:15:53Z

Another thing I've noticed is that you've based it off a fairly old commit. You'll need to rebase it onto the latest main please.

cortner · 2022-01-07T02:16:46Z

You'll see there that EuclideanVector has changed a bit and it would be good to model EuclideanMatrix on that. I hope that should require minimal work.

EquivariantMatrix -> EuclideanMatrix

MatthiasSachs · 2022-01-07T20:18:13Z

You'll see there that EuclideanVector has changed a bit and it would be good to model EuclideanMatrix on that. I hope that should require minimal work.

This is mostly done. However,

I am not sure if the implementation of coco_init(::EuclideanMatrix{CT}) is correct?
Moreover, somehow there is an error thrown which relates to the conversion between SMatrix and EuclideanMatrix.

Would you mind taking a brief look at that (I think this might be much easier for you to resolve since you know the parts of the code well), or alternatively, we could meet at the beginning of next week?

cortner · 2022-01-08T00:59:27Z

Error message:

 [1] convert(T::Type, φ::EuclideanMatrix{Float64})
   @ ACE ~/gits/temp/ms.ACE.jl/src/properties.jl:28
 [2] push!(a::Vector{EuclideanMatrix{ComplexF64}}, item::EuclideanMatrix{Float64})
   @ Base ./array.jl:994
....

It happens in this line:

push!(vals, U[irow, icol])

The computed cc U[irow, icol] is of type item::EuclideanMatrix{Float64} and get's appended to a vector with eltype EuclideanMatrix{ComplexF64}. We need a convert implementation for this. THat's easy to do, but it concerns me. Why are we creating a complex-valued CC matrix but then fill it with real-valued CCs?

MatthiasSachs · 2022-01-08T15:05:36Z

Error message:
 [1] convert(T::Type, φ::EuclideanMatrix{Float64})
   @ ACE ~/gits/temp/ms.ACE.jl/src/properties.jl:28
 [2] push!(a::Vector{EuclideanMatrix{ComplexF64}}, item::EuclideanMatrix{Float64})
   @ Base ./array.jl:994
....
It happens in this line:
push!(vals, U[irow, icol])
The computed cc U[irow, icol] is of type item::EuclideanMatrix{Float64} and get's appended to a vector with eltype EuclideanMatrix{ComplexF64}. We need a convert implementation for this. THat's easy to do, but it concerns me. Why are we creating a complex-valued CC matrix but then fill it with real-valued CCs?

I fixed other parts of the code so that now there is no conversion required! It all seems to work now.

cortner

looks mostly ok, just a few more smaller things to address please. Can you then also add the new tests to the test suite please?

The zero-correlation implementation looks roughly ok, but I added a small comment.

cortner · 2022-01-08T17:18:24Z

src/properties.jl

-real(φ::EuclideanMatrix) = EuclideanMatrix(φ.val)
-complex(φ::EuclideanMatrix) = EuclideanMatrix(φ.val)
+real(φ::EuclideanMatrix) = EuclideanMatrix(real.(φ.val))
+complex(φ::EuclideanMatrix) = EuclideanMatrix(complex(φ.val))
 complex(::Type{EuclideanMatrix{T}}) where {T} = EuclideanMatrix{complex(T)}

 +(x::SMatrix{3}, y::EuclideanMatrix) = EuclideanMatrix(x + y.val)
 Base.convert(::Type{SMatrix{3, 3, T, 9}}, φ::EuclideanMatrix) where {T} =  convert(SMatrix{3, 3, T, 9}, φ.val)


where is this needed?

not sure, I kind of copied the implementations of EuclideanVector here... Should I check, or should we introduce a new abstract type to avoid code redundancy?

ok, leave it for now.

cortner · 2022-01-08T17:19:13Z

src/properties.jl

-struct EuclideanMatrix{T} <: AbstractProperty where T<:Real
-   val::SMatrix{3, 3, Complex{T}, 9}
+struct EuclideanMatrix{T} <: AbstractProperty
+   val::SMatrix{3, 3, T, 9}
 end

 Base.show(io::IO, φ::EuclideanMatrix) = 


I'm not sure about this stype of printing; it suggests 3 vectors rather than one matrix. Would it be ok to print it similarly to the spherical matrix, but with an "e" label instead of a "y" label? (yes, I also don't like the "m" label)

Done! But I am wondering whether we should have a different label than the label of Euclidean vectors? (Since EuclideanMatrix satisfies a different equivariance property.)

I think it can be clear from the context that it is a matrix. If it were e and then three vectors then I would be confused.

cortner · 2022-01-08T17:22:58Z

src/properties.jl

+   if length(bb) == 0  # no zero-correlations allowed 
+      return false 
+   end
+   if length(bb) == 1 #MS: Not sure if this should be here


I don't remember the issue - same with EuclideanVector. If codes run ok, leave it for now but file an issue to investigate this. Do you have an idea?

The code runs. Issue is filed (issue #100)!

cortner · 2022-01-08T17:25:07Z

src/properties.jl

         : coco_zeros(phi, l, m, μ, T, A)  )


-coco_init(::EuclideanMatrix{CT}) where {CT<:Real} = [EuclideanMatrix{CT}(SMatrix{3,3,Float64,9}([1.0,0,0,0,1.0,0,0,0,1.0]))]       


I think you want

EuclideanMatrix{CT}(SMatrix{3,3,CT}([1.0,0,0,0,1.0,0,0,0,1.0]))

If I remember correctly this is the implementation for zero-correlations. So you are right that identity is the only zero-correlation there is, but you probably need to specify the right eltype.

included in issue #100

cortner · 2022-01-08T17:26:05Z

src/properties.jl

         : coco_zeros(phi, l, m, μ, T, A)  )


-coco_init(::EuclideanMatrix{CT}) where {CT<:Real} = [EuclideanMatrix{CT}(SMatrix{3,3,Float64,9}([1.0,0,0,0,1.0,0,0,0,1.0]))]       
+#coco_init(::EuclideanMatrix{CT}) where {CT<:Real} = [EuclideanMatrix(SMatrix{3,3,Complex{CT},9}([1.0,0,0,0,1.0,0,0,0,1.0]))]       
 coco_type(φ::EuclideanMatrix) = typeof(complex(φ))


so this confuses me now. In the previous commit the CCs there were produced were all real. But here you are saying they should be complex. So which is it?

This is all very confusing! But I know for sure that the coupling coefficients should be complex (same as in the implementation of EuclideanVector). Thus, my understanding is that both the coco_init function as well as the coco_zeros function should return complex valued matrices.

cortner · 2022-01-08T17:26:46Z

src/properties.jl

 coco_type(φ::EuclideanMatrix) = typeof(complex(φ))
 coco_type(::Type{EuclideanMatrix{T}}) where {T} = EuclideanMatrix{complex(T)}

 # This is slightly different from implementation in EuclideanVector!
-coco_zeros(::EuclideanMatrix, ll, mm, kk, T, A) =  EuclideanMatrix.(zeros(SMatrix{3, 3, Complex{T}, 9},9))
+coco_zeros(φ::EuclideanMatrix, ll, mm, kk, T, A) =  zeros(typeof(complex(φ)), 9)


same here, why are you enforcing complex CCs?

MatthiasSachs · 2022-01-09T16:13:35Z

sorry, just pushed a fix for a small bug that I noticed today in the morning.

- updated implementation of EuclideanVector (moved elementary rotation coefficients to dictionary and changed format) - updated dictionary for elementary rotation coefficients to inlclude l=0,1,2 All tests pass, there are now non-symmetric matrix valued basis functions, BUT the complex part is still not negligible.

-new filter added to basis selectors - utils.SymmetricBond_basis extended to include symmetry option for bond sign flip - renamed "test_basisselectors.jl"

MatthiasSachs · 2022-02-02T18:10:13Z

@cortner I think this PR is ready to merge. Do you approve?

cortner · 2022-02-02T19:45:48Z

I will take a look. We will also need to merge in the latest main branch, or rebase.

…ices

MatthiasSachs · 2022-02-03T23:05:24Z

@cortner: I just merged the upstream/main into this branch.

MatthiasSachs · 2022-02-03T23:22:08Z

Argghh... it's not obvious to me how to fix the newly appearing error after merging with upstream/main!

cortner · 2022-02-04T21:14:18Z

I wil try - no promises about time-scale

cortner · 2022-02-04T21:50:12Z

what error? I only get a warning?

MatthiasSachs · 2022-02-04T21:54:26Z

This one:

ERROR: LoadError: MethodError: namedtuple(::Tuple{}) is ambiguous. Candidates:
  namedtuple(names::Tuple{Vararg{Symbol, N}}) where N in NamedTupleTools at /Users/msachs2/.julia/packages/NamedTupleTools/GZxRn/src/NamedTupleTools.jl:196
  namedtuple(names::Tuple{Vararg{String, N}}) where N in NamedTupleTools at /Users/msachs2/.julia/packages/NamedTupleTools/GZxRn/src/NamedTupleTools.jl:198
Possible fix, define
  namedtuple(::Tuple{})
Stacktrace:
  [1] select(nt::NamedTuple{(:be, :n, :l, :m), Tuple{Symbol, Int64, Int64, Int64}}, ks::Tuple{})
    @ NamedTupleTools ~/.julia/packages/NamedTupleTools/GZxRn/src/NamedTupleTools.jl:270
  [2] _broadcast_getindex_evalf
    @ ./broadcast.jl:648 [inlined]
  [3] _broadcast_getindex
    @ ./broadcast.jl:621 [inlined]
  [4] getindex
    @ ./broadcast.jl:575 [inlined]
  [5] copy
    @ ./broadcast.jl:922 [inlined]
  [6] materialize
    @ ./broadcast.jl:883 [inlined]
  [7] _sparsify_component!(basis1p::ACE.EllipsoidBondEnvelope{Float64}, keep::Vector{NamedTuple{(:be, :n, :l, :m), Tuple{Symbol, Int64, Int64, Int64}}})
    @ ACE ~/.julia/dev/ACE/src/product_1pbasis.jl:297
  [8] sparsify!(basis1p::ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, keep::Vector{NamedTuple})
    @ ACE ~/.julia/dev/ACE/src/product_1pbasis.jl:277
  [9] clean_1pbasis!(basis::PIBasis{ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, typeof(real), Float64, ComplexF64})
    @ ACE ~/.julia/dev/ACE/src/pibasis.jl:221
 [10] clean_pibasis!(basis::ACE.SymmetricBasis{PIBasis{ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, typeof(real), Float64, ComplexF64}, ACE.Invariant{Float64}, ACE.O3{:l, :m}, typeof(real), ACE.Invariant{Float64}}; atol::Float64)
    @ ACE ~/.julia/dev/ACE/src/symmbasis.jl:250
 [11] clean_pibasis!(basis::ACE.SymmetricBasis{PIBasis{ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, typeof(real), Float64, ComplexF64}, ACE.Invariant{Float64}, ACE.O3{:l, :m}, typeof(real), ACE.Invariant{Float64}})
    @ ACE ~/.julia/dev/ACE/src/symmbasis.jl:243
 [12] ACE.SymmetricBasis(φ::ACE.Invariant{Float64}, symgrp::ACE.O3{:l, :m}, pibasis::PIBasis{ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, typeof(real), Float64, ComplexF64}, _real::Function)
    @ ACE ~/.julia/dev/ACE/src/symmbasis.jl:176
 [13] #SymmetricBasis#313
    @ ~/.julia/dev/ACE/src/symmbasis.jl:103 [inlined]
 [14] ACE.SymmetricBasis(φ::ACE.Invariant{Float64}, basis1p::ACE.Product1pBasis{4, Tuple{Categorical1pBasis{:be, :be, 2, Symbol}, ACE.Rn1pBasis{Float64, PolyTransform{Int64, Float64}, ACE.OrthPolys.OrthPolyBasis{Float64}, :rr, :n, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, Float64}}}}}, ACE.Ylm1pBasis{Float64, :rr, :l, :m, ACE.DState{NamedTuple{(:rr,), Tuple{SVector{3, ComplexF64}}}}}, ACE.EllipsoidBondEnvelope{Float64}}, ComplexF64}, symgrp::ACE.O3{:l, :m}, Bsel::ACE.CategorySparseBasis; isreal::Bool, kwargs::Base.Iterators.Pairs{Symbol, ACE.Utils.var"#5#7", Tuple{Symbol}, NamedTuple{(:filterfun,), Tuple{ACE.Utils.var"#5#7"}}})
    @ ACE ~/.julia/dev/ACE/src/symmbasis.jl:96
 [15] #SymmetricBasis#310
    @ ~/.julia/dev/ACE/src/symmbasis.jl:85 [inlined]
 [16] SymmetricBond_basis(ϕ::ACE.Invariant{Float64}, env::ACE.EllipsoidBondEnvelope{Float64}, Bsel::SparseBasis; RnYlm::Nothing, bondsymmetry::Nothing, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ACE.Utils ~/.julia/dev/ACE/src/utils.jl:97
 [17] SymmetricBond_basis(ϕ::ACE.Invariant{Float64}, env::ACE.EllipsoidBondEnvelope{Float64}, Bsel::SparseBasis)
    @ ACE.Utils ~/.julia/dev/ACE/src/utils.jl:78
 [18] top-level scope
    @ ~/.julia/dev/ACE/test/test_bondbasisselectors.jl:38
in expression starting at /Users/msachs2/.julia/dev/ACE/test/test_bondbasisselectors.jl:36

This is the error I get on my own computer. But I presume that the error encountered in 4b82aec is the same.

cortner · 2022-02-04T22:20:00Z

ok, I can reproduce

cortner · 2022-02-04T22:26:37Z

Ok - the problem is here:

env = ACE.EllipsoidBondEnvelope(r0cut, rcut, zcut;floppy=false, λ= .5)
property = Invariant() ... 
basis = ACE.Utils.SymmetricBond_basis(property, env, Bsel; )

you are trying to create a symmetric basis but it isn't a basis, so the basis spec contains only empty tuples. That's a weird edge case and I'm not even sure I want to allow this... Is this every going to occur outside of a test?

cortner · 2022-02-04T22:27:38Z

See here : once we sparsity (in this case just "clean up") ...

function _sparsify_component!(basis1p, keep)
   # get rid of all info we don't need 
   syms = symbols(basis1p)
   keep1 = unique( select.(keep, Ref(syms)) )
   ....

I can't "select" a subset of the basis specs because there is no spec. It is at odds with the entire design of the basis framework.

MatthiasSachs · 2022-02-05T20:54:56Z

Ok - the problem is here:
env = ACE.EllipsoidBondEnvelope(r0cut, rcut, zcut;floppy=false, λ= .5)
property = Invariant() ... 
basis = ACE.Utils.SymmetricBond_basis(property, env, Bsel; )
you are trying to create a symmetric basis but it isn't a basis, so the basis spec contains only empty tuples. That's a weird edge case and I'm not even sure I want to allow this... Is this every going to occur outside of a test?

I can see in the error message that the spec is empty but I don't yet understand why this is the case? I create a special basis selector within the call of SymmetricBond_basis and then call a constructor for SymmetricBasis using this basis selector. Why is what I am instantiating not a basis? Is the bond envelope env and the fact that we don't associate a symbol with that causing problems here? So far this is the standard way for me to create bond bases...

cortner · 2022-02-05T22:05:12Z

It’s just a multiplier by itself - it isn’t multiplying a basis

cortner · 2022-02-15T05:31:02Z

Any thoughts on this? Maybe there is a way to just change the tests? Or do you need me to figure out a way around this? I'M just a bit overwhelmed right now - sorry for the delay.

Matthias Sachs and others added 2 commits October 29, 2021 11:25

tests passing for tol = 1e-5

696e79e

more accurate rotation coefficients added

e6ba036

cortner and others added 4 commits January 6, 2022 18:30

fixed fio for EquivariantMatrix

884ebf9

renaming

8dc3526

EquivariantMatrix -> EuclideanMatrix

Merge branch 'main' into ms/equivarinat_matrices

75d5b6a

rebased branch and adapted implementation

d1fa412

conversion issue fixed. All tests pass!

f81c451

added test_EuclideanMatrix.jl to test set

c085a44

cortner requested changes Jan 8, 2022

View reviewed changes

MatthiasSachs added 3 commits January 8, 2022 17:28

derivative test removed

f445e56

show modified and test for complex version of EuclideanMatrix added

9dcc9bb

small fix

f6ea6e8

additional tests added (some don't pass!)

c7cdea1

MatthiasSachs mentioned this pull request Jan 11, 2022

Weird issues with EuclideanMatrix property #104

Open

cortner and others added 6 commits January 11, 2022 15:56

small changes to show and testing

fa885f2

several changes

28ac1fb

-new filter added to basis selectors - utils.SymmetricBond_basis extended to include symmetry option for bond sign flip - renamed "test_basisselectors.jl"

Merge branch 'main' into ms/equivarinat_matrices

b035504

small change to tests

48bca83

clean-up in tests for bond basis selectors

0991821

update

a06a303

MatthiasSachs added 2 commits February 2, 2022 18:33

Documentation of filter function EvenL updated

e856741

edit fixed

e02dd36

Merge remote-tracking branch 'upstream/main' into ms/equivarinat_matr…

4b82aec

…ices

slim printing of test

6793aae

cortner and others added 3 commits March 3, 2022 10:31

fix sparsification for bses without symbol

1501aca

added mixed license text

8e178e0

license stuff

298c8c3

cortner merged commit f148042 into ACEsuit:main Mar 3, 2022

		: coco_zeros(phi, l, m, μ, T, A) )


		coco_init(::EuclideanMatrix{CT}) where {CT<:Real} = [EuclideanMatrix{CT}(SMatrix{3,3,Float64,9}([1.0,0,0,0,1.0,0,0,0,1.0]))]

Ms/equivarinat matrices #99

Ms/equivarinat matrices #99

Conversation

MatthiasSachs commented Jan 5, 2022

MatthiasSachs commented Jan 5, 2022

cortner commented Jan 5, 2022 • edited Loading

cortner commented Jan 7, 2022

cortner commented Jan 7, 2022

cortner commented Jan 7, 2022

cortner commented Jan 7, 2022

MatthiasSachs commented Jan 7, 2022

cortner commented Jan 8, 2022

MatthiasSachs commented Jan 8, 2022

cortner left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MatthiasSachs Jan 8, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MatthiasSachs Jan 8, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MatthiasSachs commented Jan 9, 2022 • edited Loading

MatthiasSachs commented Feb 2, 2022

cortner commented Feb 2, 2022

MatthiasSachs commented Feb 3, 2022

MatthiasSachs commented Feb 3, 2022 • edited Loading

cortner commented Feb 4, 2022

cortner commented Feb 4, 2022

MatthiasSachs commented Feb 4, 2022 • edited Loading

cortner commented Feb 4, 2022

cortner commented Feb 4, 2022

cortner commented Feb 4, 2022

MatthiasSachs commented Feb 5, 2022

cortner commented Feb 5, 2022

cortner commented Feb 15, 2022 • edited Loading

cortner commented Jan 5, 2022 •

edited

Loading

MatthiasSachs Jan 8, 2022 •

edited

Loading

MatthiasSachs Jan 8, 2022 •

edited

Loading

MatthiasSachs commented Jan 9, 2022 •

edited

Loading

MatthiasSachs commented Feb 3, 2022 •

edited

Loading

MatthiasSachs commented Feb 4, 2022 •

edited

Loading

cortner commented Feb 15, 2022 •

edited

Loading