Coverage for dynasor/modes/tools.py: 57%

1import numbers

2from fractions import Fraction

3from typing import Optional, Union

4import numba

5import numpy as np

6from numpy.typing import NDArray

9# Linalg

10def trace(A):

11 return A[0, 0] + A[1, 1] + A[2, 2]

14def det(A):

15 if len(A) == 2:

16 return A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0]

17 d = 0

18 for i, B in enumerate(A[0]):

19 minor = np.hstack([A[1:, :i], A[1:, i+1:]])

20 d += (-1)**i * B * det(minor)

21 return d

24def inv(

25 A: Union[Fraction, NDArray[float]],

26 as_fraction: Optional[bool] = True,

27) -> Union[Fraction, NDArray[float]]:

28 """

29 Inverts the A matrix.

31 Parameters

32 ----------

33 A

34 Array to be inverted.

35 as_fraction

36 Boolean which determines if inverted matrix is returned

37 as :class:`Fraction` matrix or real matrix.

38 """

40 detx2 = det(A) * 2

42 rest = (trace(A)**2 - trace(A @ A)) * np.diag([1, 1, 1]) - 2 * A * trace(A) + 2 * A @ A

44 if detx2 < 0:

45 detx2 = -detx2

46 rest = -rest

48 A_inv = np.array([Fraction(n, detx2) for n in rest.flat]).reshape(3, 3)

50 assert np.all(A @ A_inv == np.eye(len(A)))

52 # return either real matrix or as fraction matrix

53 A_inv_real = np.reshape([float(x) for x in A_inv.flat], A_inv.shape)

54 assert np.allclose(A_inv_real, np.linalg.inv(A))

56 if as_fraction: 56 ↛ 59line 56 didn't jump to line 59 because the condition on line 56 was always true

57 return A_inv

58 else:

59 return A_inv_real

62def symmetrize_eigenvectors(

63 eigenvectors: NDArray[float],

64 cell: Optional[NDArray[float]] = None,

65 max_iter: Optional[int] = 1000,

66 method: Optional[str] = 'varimax',

67) -> NDArray[float]:

68 """Takes a set of vectors and tries to make them nice

70 Parameters

71 ----------

72 eigenvectors

73 If there are n degenerate eigenvectors and m atoms in the basis the

74 array should be `(n,m,3)`.

75 cell

76 If default `None` nothing is done to the cartesian directions but a cell

77 can be provided so the directions are in scaled coordinates instead.

78 max_iter

79 The number of iterations in the symmetrization procedure.

80 method

81 Can be `'varimax'` or `'quartimax'` or a parameter between `0: qartimax` and

82 `1: varimax`. Depending on the choice one obtains, e.g., Equamax, Parsimax, etc.

84 """

86 if cell is None: 86 ↛ 90line 86 didn't jump to line 90 because the condition on line 86 was always true

87 cell = np.eye(3)

89 # s = band, i = basis, a = axis

90 eigenvectors = np.einsum('sia,ab->sib', eigenvectors, np.linalg.inv(cell))

92 components = eigenvectors.reshape(len(eigenvectors), -1).T

94 rotation_matrix = factor_analysis(components, iterations=max_iter, method=method)

96 new_eigenvectors = np.dot(components, rotation_matrix).T

98 new_eigenvectors = new_eigenvectors.reshape(len(new_eigenvectors), -1, 3)

100 new_eigenvectors = np.einsum('sia,ab->sib', new_eigenvectors, cell)

101

102 return new_eigenvectors

103

104

105def factor_analysis(

106 L: NDArray[float],

107 iterations: Optional[int] = 1000,

108 method: Optional[str] = 'varimax',

109) -> NDArray[float]:

110 """Performs factor analysis on `L` finding rotation matrix `R` such that `L @ R = L'` is simple.

111

112 In the future consider using the scikit learn methods directly but beware

113 the changes need to accomodate complex numbers.

114

115 References:

116 * *Sparse Modeling of Landmark and Texture Variability using the Orthomax Criterion*

117 Mikkel B. Stegmann, Karl Sjöstrand, Rasmus Larsen

118 http://www2.imm.dtu.dk/pubdb/edoc/imm4041.pdf

119

120 * http://www.cs.ucl.ac.uk/staff/d.barber/brml

121

122 * https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.FactorAnalysis.html # noqa

123

124 * https://stats.stackexchange.com/questions/185216/factor-rotation-methods-varimax-quartimax-oblimin-etc-what-do-the-names # noqa

125

126 """

127

128 nrow, ncol = L.shape

129 R = np.eye(ncol)

130

131 if method == 'varimax': 131 ↛ 133line 131 didn't jump to line 133 because the condition on line 131 was always true

132 gamma = 1

133 elif method == 'quartimax':

134 gamma = 0

135 else:

136 gamma = method

137

138 for _ in range(iterations):

139 LR = np.dot(L, R)

140 grad = LR * (np.abs(LR)**2 - gamma * np.mean(np.abs(LR)**2, axis=0))

141 G = L.T.conj() @ grad

142 u, s, vh = np.linalg.svd(G)

143 R = u @ vh

144

145 return R

146

147

148def group_eigvals(

149 vals: NDArray[float],

150 tol: Optional[float] = 1e-6,

151) -> list[list[int]]:

152 assert sorted(vals) == list(vals), vals

153

154 groups = [[0]]

155 for i in range(1, len(vals)):

156 if np.abs(vals[i] - vals[i-1]) < tol:

157 groups[-1].append(i)

158 else:

159 groups.append([i])

160 return groups

161

162

163# misc

164def as_fraction(not_fraction):

165 if isinstance(not_fraction, numbers.Number):

166 return Fraction(not_fraction)

167

168 if isinstance(not_fraction, np.ndarray):

169 arr = np.array([Fraction(n) for n in not_fraction.flat])

170 return arr.reshape(not_fraction.shape)

171

172 if isinstance(not_fraction, tuple):

173 arr = (Fraction(n) for n in not_fraction)

174 return arr

175

176 if isinstance(not_fraction, list):

177 arr = [Fraction(n) for n in not_fraction]

178 return arr

179

180

181@numba.njit

182def get_dynamical_matrix(fc, offsets, indices, q):

183

184 n = indices.max() + 1

185 N = len(fc)

186 D = np.zeros(shape=(n, n, 3, 3), dtype=np.complex128)

187 for ia in range(n):

188 for a in range(N):

189 if ia != indices[a]:

190 continue

191 na = offsets[a]

192 for b in range(N):

193 ib = indices[b]

194 nb = offsets[b]

195 dn = (nb - na).astype(np.float64)

196 D[ia, ib] += fc[a, b] * np.exp(2j*np.pi * np.dot(dn, q))

197 break

198 return D

199

200

201# For debug, this is a slower but perhaps more accurate variant, if the fc

202# obeys translational invariance this should give the same result

203# @numba.njit

204# def get_dynamical_matrix_full(fc, offsets, indices, q):

205#

206# n = indices.max() + 1

207# N = len(fc)

208# D = np.zeros(shape=(n, n, 3, 3), dtype=np.complex128)

209# for I in range(N):

210# i = indices[I]

211# m = offsets[I]

212# for J in range(N):

213# j = indices[J]

214# n = offsets[J]

215#

216# off = (n - m).astype(np.float64)

217#

218# phase = np.exp(1j * 2*np.pi * np.dot(off, q))

219#

220# D[i, j] += fc[I, J] * phase

221#

222# D /= (N / D.shape[0])

223#

224# return D

225

226

227def make_table(M):

228 rows = []

229 for r in M:

230 rows.append(''.join(f'{e:<20.2f}' for e in r))

231 return '\n'.join(rows)

Coverage for dynasor / modes / tools.py: 57%

96 statements