-
Notifications
You must be signed in to change notification settings - Fork 185
/
K_symmetrised.m
114 lines (92 loc) · 3 KB
/
K_symmetrised.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
function w = K_symmetrised(psi,q1,q2,CS,SS,varargin)
% evaluate kernel modulo symmetries
%
% Input
% psi - @SO3Kernel
% q1, q2 - @quaternion(s)
% CS, SS - crystal , specimen @symmetry
%
% Options
% exact - Use complete kernel functions and compute full distance matrix
% epsilon - cut kernel functions outside of [-2epsilon,2epsilon]
%
% Description
% K(q1,q2) = Sum(S) Sum(l) A_l Tr T_l(s1^-1 q1 s2)
%
% See also
% SO3FunRBF/eval SO3Fun/interpolate
% only the pur rotational part is of interest
% TODO
qCS = unique(quaternion(CS));
qSS = unique(quaternion(SS));
if check_option(varargin,'exact')
epsilon = pi;
else
epsilon = min(pi,get_option(varargin,'epsilon',psi.halfwidth*3.5));
end
% how to use sparse matrix representation
if isa(q1,'SO3Grid')
lg1 = length(q1);
else
lg1 = -length(q1);
end
if isa(q2,'SO3Grid')
lg2 = length(q2);
else
lg2 = -length(q2);
end
% TODO: Better use the diameter here fundamentalRegion(CS,SS).maxAngle
% For triclinic symmetries the right hand side should be pi and >=
% Note that the following condition occurs in SO3Fun.interpolate
if epsilon>2*pi/CS.Laue.multiplicityZ % full matrices
q1 = quaternion(q1);
q2 = quaternion(q2);
w = zeros(length(q1),length(q2));
for iks = 1:length(qCS)
for ips = 1:length(qSS) % for all symmetries
sg = qSS(ips) * q1 * qCS(iks); % rotate g1
omega = abs(dot_outer(sg,q2)); % calculate full distance matrix
w = w + psi.eval(omega);
end
end
elseif (lg1>0 || lg2>0) && ~check_option(varargin,'old')
w = sparse(abs(lg1),abs(lg2));
% sum over specimen symmetry
if (lg1 >= lg2) % first argument is SO3Grid
for issq = 1:length(qSS)
d = abs(dot_outer(q1,qSS(issq)*quaternion(q2),'epsilon',epsilon,...
'nospecimensymmetry'));
w = w + spfun(@psi.eval,d);
end
else % second argument is SO3Grid
for issq = 1:length(qSS)
d = abs(dot_outer(q2,qSS(issq)*quaternion(q1),'epsilon',epsilon,...
'nospecimensymmetry'));
w = w + spfun(@psi.eval,d.');
end
end
else
q1 = quaternion(q1);
q2 = quaternion(q2);
w = sparse(length(q1),length(q2));
for iks = 1:length(qCS)
for ips = 1:length(qSS) % for all symmetries
if abs(lg1) > abs(lg2)
sg = qSS(ips) * q2 * qCS(iks); % rotate g1
omega = abs(dot_outer(q1,sg)); % calculate full distance matrix
else
sg = qSS(ips) * q1 * qCS(iks); % rotate g1
omega = abs(dot_outer(sg,q2)); % calculate full distance matrix
end
% z = find(omega>cos(epsilon));
% if length(z) > length(omega)/length(CS)/10, w = full(w); end
% w(z) = w(z) + kk.K(omega(z));
% if length(z) > numel(omega)/length(CS)/10, w = full(w); end
[y,x] = find(omega>cos(epsilon));
dummy = sparse(y,x,psi.eval(omega(sub2ind(size(w),y,x))),length(q1),length(q2));
w = w + dummy;
end
end
end
%nnz(w)
w = w / length(qCS) / length(qSS);