15 KiB
Testing Svd Testing SVD
- Generate a block random distribution of positions
- Is the SVD decomposition of distances really useful ?
- Generating a box of electrons centered around nucleii
- Calculating SVD of the distance vector
Generate a block random distribution of positions
Here's a small function to generate random points inside a 3D box of given length.
def get_sphere_distribution(n, dmin, Ls, maxiter=1e4, allow_wall=True):
"""Get random points in a box with given dimensions and minimum separation.
Parameters:
- n: number of points
- dmin: minimum distance
- Ls: dimensions of box, shape (3,) array
- maxiter: maximum number of iterations.
- allow_wall: whether to allow points on wall;
(if False: points need to keep distance dmin/2 from the walls.)
Return:
- ps: array (n, 3) of point positions,
with 0 <= ps[:, i] < Ls[i]
- n_iter: number of iterations
- dratio: average nearest-neighbor distance, divided by dmin.
Note: with a fill density (sphere volume divided by box volume) above about
0.53, it takes very long. (Random close-packed spheres have a fill density
of 0.64).
Author: Han-Kwang Nienhuys (2020)
Copying: BSD, GPL, LGPL, CC-BY, CC-BY-SA
See Stackoverflow: https://stackoverflow.com/a/62895898/6228891
"""
Ls = np.array(Ls).reshape(3)
if not allow_wall:
Ls -= dmin
# filling factor; 0.64 is for random close-packed spheres
# This is an estimate because close packing is complicated near the walls.
# It doesn't work well for small L/dmin ratios.
sphere_vol = np.pi/6*dmin**3
box_vol = np.prod(Ls + 0.5*dmin)
fill_dens = n*sphere_vol/box_vol
if fill_dens > 0.64:
msg = f'Too many to fit in the volume, density {fill_dens:.3g}>0.64'
raise ValueError(msg)
# initial try
ps = np.random.uniform(size=(n, 3)) * Ls
# distance-squared matrix (diagonal is self-distance, don't count)
dsq = ((ps - ps.reshape(n, 1, 3))**2).sum(axis=2)
dsq[np.arange(n), np.arange(n)] = np.infty
for iter_no in range(int(maxiter)):
# find points that have too close neighbors
close_counts = np.sum(dsq < dmin**2, axis=1) # shape (n,)
n_close = np.count_nonzero(close_counts)
if n_close == 0:
break
# Move the one with the largest number of too-close neighbors
imv = np.argmax(close_counts)
# new positions
newp = np.random.uniform(size=3)*Ls
ps[imv]= newp
# update distance matrix
new_dsq_row = ((ps - newp.reshape(1, 3))**2).sum(axis=-1)
dsq[imv, :] = dsq[:, imv] = new_dsq_row
dsq[imv, imv] = np.inf
else:
raise RuntimeError(f'Failed after {iter_no+1} iterations.')
if not allow_wall:
ps += dmin/2
dratio = (np.sqrt(dsq.min(axis=1))/dmin).mean()
return ps, iter_no+1, dratio
A function that generates a block random set of points
<<generate3Dpoints>>
def generateBlockRandomPointsAtShftApart(n,L1,dmin,shift):
"""
Parameters:
- n: number of points
- L1: dimensions of box, shape (3,) array
- dmin: minimum intra block distance
- shift: inter block distance
Return:
- r: array (n, 3) of point positions,
"""
Ls = np.array([L1,L1,L1]) # lengths of the box
# Quadrant +,+
dist11,a,b = get_sphere_distribution(n,dmin,Ls)
dist11 += np.array([shift,shift,0.0])
# Quadrant -,+
dist21,a,b = get_sphere_distribution(n,dmin,Ls)
dist21 += np.array([-shift,shift,0.0])
# Quadrant +,-
dist12,a,b = get_sphere_distribution(n,dmin,Ls)
dist12 += np.array([shift,-shift,0.0])
# Quadrant -,-
dist22,a,b = get_sphere_distribution(n,dmin,Ls)
dist22 += np.array([-shift,-shift,0.0])
xs = np.concatenate((dist11[:,0],dist12[:,0],dist21[:,0],dist22[:,0]))
ys = np.concatenate((dist11[:,1],dist12[:,1],dist21[:,1],dist22[:,1]))
zs = np.concatenate((dist11[:,2],dist12[:,2],dist21[:,2],dist22[:,2]))
return((np.array((xs,ys,zs))).T)
None
Is the SVD decomposition of distances really useful ?
The problem is the following:
Given a metric F(r1,r2) finding the best representatin of F(r1,r2)
Slater metric
The metric can be a scaled distance such as
\[ F(r_1,r_2) = \exp(-\kappa(|r_1 - r_2|)) \]
Where,\(\kappa\) is the rate of the breadth of the slater.
Gaussian metric
The metric can be a scaled distance such as
\[ F(r_1,r_2) = \exp(-\kappa(|r_1 - r_2|^2)) \]
Where,\(\kappa\) is the rate of the breadth of the gaussian.
Generating a box of electrons centered around nucleii
Generate randomly distributed nucleii
The nucleii will serve as anchors for the screening of distances for calculating the screened \(\tilde{F}(r_1,r_2)\) metric.
Calculating SVD of the distance vector
Slater metric
Calculate the matrix of the \(F(r_1,r_2)\) metric i.e. the slater metric.
import numpy
a = numpy.array([[1,2,3],[4,5,6],[7,8,9]])
b = numpy.array([[11,12,13],[14,15,16],[17,18,19]])
print(list(zip(a,b))[0][1])
print(numpy.square(a[:,0]))
def stepExp(a):
def myexp(x):
if numpy.abs(x) > 1e+0:
return numpy.zeros_like(x)
else:
return numpy.exp(x)
res = numpy.array([[myexp(x) for x in y] for y in a])
return(res)
print(numpy.exp(a))
print(stepExp(a))
[11 12 13] [ 1 16 49] [[2.71828183e+00 7.38905610e+00 2.00855369e+01] [5.45981500e+01 1.48413159e+02 4.03428793e+02] [1.09663316e+03 2.98095799e+03 8.10308393e+03]] [[2.71828183 0. 0. ] [0. 0. 0. ] [0. 0. 0. ]]
Gaussian metric
Calculate the matrix of the \(FG(r_1,r_2)\) metric i.e. the gaussian metric.
Palying around
Calculate the matrix of the \(FG(r_1,r_2)\) metric i.e. the gaussian metric.
Testing SVD for custom matrices
import numpy
import matplotlib.pyplot as plt
a = numpy.array([[0,100,200],[100,0,200],[100,200,0]])
b = numpy.exp(-a)
print("Matrix A")
print(a)
print("Matrix Exp(A)")
print(numpy.around(b,10))
u,d,vt = numpy.linalg.svd(a)
d = d/numpy.linalg.norm(d)
print("Singular values of A")
print(numpy.around(d,3))
print("Singular vectors of A")
print(numpy.around(u,3))
u,d,vt = numpy.linalg.svd(b)
d = d/numpy.linalg.norm(d)
print("Singular values of Exp(A)")
print(numpy.around(d,3))
print("Singular vectors of Exp(A)")
print(numpy.around(u,3))
Matrix A [[ 0 100 200] [100 0 200] [100 200 0]] Matrix Exp(A) [[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]] Singular values of A [0.813 0.53 0.24 ] Singular vectors of A [[-0.67 0.142 0.728] [-0.626 0.42 -0.657] [-0.399 -0.896 -0.193]] Singular values of Exp(A) [0.577 0.577 0.577] Singular vectors of Exp(A) [[-1. 0. -0. ] [-0. -0.894 -0.447] [-0. -0.447 0.894]]