basis functions are not consistent with different MaxDeg

[hjchen1983]

• julia Version 1.7.3
• ACE Version v0.12.43

  
  using ACE, LinearAlgebra
  using ACE: scal1pbasis, λ, PIBasis, SymmetricBasis, LinearACEModel,
         evaluate, grad_config, grad_params
  const ↑ = Symbol('↑')
  const ↓ = Symbol('↓')
  Nels = 1
  trans = λ("r -> atan(r)")
  x0 = -3.0
  x1 =  3.0
  rand_el() = State(x = rand()*(x1-x0) + x0, σ = rand([↑, ↓]))
  cfg = [ rand_el() for i in 1:Nels ]

maxdeg = 8 P = scal1pbasis(:x, :n, maxdeg, trans, x1, x0; pin = 1, pcut = 1) P.meta["x0"] = x0; P.meta["x1"] = x1 S = Categorical1pBasis([↑, ↓], :σ, :s, "S") B1p = ACE.Product1pBasis( (P, S) ) Bsel = SimpleSparseBasis(Nels, maxdeg) basis = PIBasis(B1p, Bsel) symbasis = SymmetricBasis(ACE.Invariant(), ACE.NoSym(), basis) s = max.(1, ACE.scaling(symbasis, 2)) c = randn(length(symbasis)) ./ s ace = LinearACEModel(symbasis, c) u0 = ACE.evaluate(ace.basis, cfg)

maxdeg = 6 P = scal1pbasis(:x, :n, maxdeg, trans, x1, x0; pin = 1, pcut = 1) P.meta["x0"] = x0; P.meta["x1"] = x1 S = Categorical1pBasis([↑, ↓], :σ, :s, "S") B1p = ACE.Product1pBasis( (P, S) ) Bsel = SimpleSparseBasis(Nels, maxdeg) basis = PIBasis(B1p, Bsel) symbasis = SymmetricBasis(ACE.Invariant(), ACE.NoSym(), basis) s = max.(1, ACE.scaling(symbasis, 2)) c = randn(length(symbasis)) ./ s ace = LinearACEModel(symbasis, c) u1 = ACE.evaluate(ace.basis, cfg)

u0 = [u0[i].val for i = 1:length(u0)] u1 = [u1[i].val for i = 1:length(u1)]

println("error of basis functions ", norm(u0[1:length(u1)] - u1));

[cortner]

Thank you. There are two possibilities: (1) the bases are genuinely inconsistent. I would consider this a bug. (2) the larger basis doesn’t respect the ordering of the smaller basis i.e. the first few basis functions just aren’t the same.

Before doing anything else can we find out which of the two it is? Can you use ACE.get_spec to look at the basis specifications for the two bases and compare values for the same basis functions?

If it turns out that (1) is the issue then I’m afraid you need to wait until I have a working laptop back. If it is (2) then you can use the basis specfications to construct a mapping (sparse matrix with just 0, 1 entries) to embed the short basis into the long one. In fact, this would be a great feature to have in ACE.jl … and I’ll think about how to do this in general.

[hjchen1983]

We think it is the first. Here is the script we use to look at the basis specifications.

using ACE, LinearAlgebra
using ACE: scal1pbasis, λ, PIBasis, SymmetricBasis, LinearACEModel,
           evaluate, grad_config, grad_params
const ↑ = Symbol('↑')
const ↓ = Symbol('↓')
Nels = 1
trans = λ("r -> atan(r)")
x0 = -3.0
x1 =  3.0
rand_el() = State(x = rand()*(x1-x0) + x0, σ = rand([↑, ↓]))
cfg = [ rand_el() for i in 1:Nels ]

maxdeg = 4
P = scal1pbasis(:x, :n, maxdeg, trans, x1, x0; pin = 1, pcut = 1)
P.meta["x0"] = x0; P.meta["x1"] = x1
S = Categorical1pBasis([↑, ↓], :σ, :s, "S")
B1p = ACE.Product1pBasis( (P, S) )
Bsel = SimpleSparseBasis(Nels, maxdeg)
basis = PIBasis(B1p, Bsel)
symbasis = SymmetricBasis(ACE.Invariant(), ACE.NoSym(), basis)
s = max.(1, ACE.scaling(symbasis, 2))
c = randn(length(symbasis)) ./ s
u0 = LinearACEModel(symbasis, c)
u0_val = ACE.evaluate(u0.basis,cfg)

maxdeg = 3
P = scal1pbasis(:x, :n, maxdeg, trans, x1, x0; pin = 1, pcut = 1)
P.meta["x0"] = x0; P.meta["x1"] = x1
S = Categorical1pBasis([↑, ↓], :σ, :s, "S")
B1p = ACE.Product1pBasis( (P, S) )
Bsel = SimpleSparseBasis(Nels, maxdeg)
basis = PIBasis(B1p, Bsel)
symbasis = SymmetricBasis(ACE.Invariant(), ACE.NoSym(), basis)
s = max.(1, ACE.scaling(symbasis, 2))
c = randn(length(symbasis)) ./ s
u1 = LinearACEModel(symbasis, c)
u1_val = ACE.evaluate(u1.basis,cfg)

# get_spec
u0_spec = ACE.get_spec(u0.basis)
u1_spec = ACE.get_spec(u1.basis)

# embedding
u0_iA = u0.basis.pibasis.spec.iAA2iA
u1_iA = u1.basis.pibasis.spec.iAA2iA

function findcol(M, col)                
    @inbounds @views for c in axes(M, 1)
        M[c,:] == col && return c       
    end                                 
    return nothing                      
end   

ind = [findcol(u0_iA,u1_iA[i,:]) for i in axes(u1_iA,1)];
ind = ind[ind .!=nothing]
'''

[cortner]

That doesn't look right. If the ordering in the A Badis is different then this won't help you.... though in the other hand you only have a single 1p basis component? So maybe this is equivalent.

Could you try to compare the output of the two P bases?

[cortner]

So I've narrowed this down to the constructioin of the discrete Jacobi polynomials:

using ACE, LinearAlgebra
using ACE: evaluate 
using ACE.OrthPolys: discrete_jacobi
J8 = discrete_jacobi(8)
J6 = discrete_jacobi(6)
x = rand() - 0.5
p8 = evaluate(J8, x)
p6 = evaluate(J6, x)
@show norm(p8[1:6] - p6)

Hopefully this is now sufficiently "local" that I can quickly fix it. But it is a scary little bug. Thank you for pointing it out!!

[cortner]

Hm - I found it. And I'm not sure anymore that it's a bug :). But it is definitely a problem. When you are constructing the basis you implicitly call this:

function discrete_jacobi(N; pcut=0, xcut=1.0, pin=0, xin=-1.0, Nquad = 3 * N, 
                            trans = identity)
...

Now the problem is the number of quadrature points that specify the orthogonality relation depends on the polynomial degree. And this is where it comes from.

My proposal is this: I change this to Nquad = max(300, 3 * N). At least this should solve your issue for now. And I'll file a new issue to remind myself that I need to return to this,

[cortner]

🎉

error of basis functions 0.0

[cortner]

this is now included in ACE.jl@0.12.44 - please change your version requirement to that lower bound.