-
Notifications
You must be signed in to change notification settings - Fork 1
/
do_schur.m
155 lines (147 loc) · 6.52 KB
/
do_schur.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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
function [ X, RR ] = do_schur( P, sd, fileid, b )
% This file is part of GPCCA.
%
% Copyright (c) 2018, 2017 Bernhard Reuter
%
% If you use this code or parts of it, cite the following reference:
%
% Reuter, B., Weber, M., Fackeldey, K., Röblitz, S., & Garcia, M. E. (2018). Generalized
% Markov State Modeling Method for Nonequilibrium Biomolecular Dynamics: Exemplified on
% Amyloid β Conformational Dynamics Driven by an Oscillating Electric Field. Journal of
% Chemical Theory and Computation, 14(7), 3579–3594. https://doi.org/10.1021/acs.jctc.8b00079
%
% GPCCA is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published
% by the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU Lesser General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------
% Perform Schur decomposition and Gram-Schmidt orthogonalization.
%
% [ X, RR ] = do_schur( P, sd, fileid )
%
% Input:
% P stochastic (fine) transition (N,N)-matrix
% sd "initial" distribution
% fileid ID-string for the naming of output files
% b if b < 0 then -b blocks will be sorted,
% if b > 0 then b or b+1 eigenvalues will be sorted,
% depending on the sizes of the blocks,
% if b = 0 then the whole Schur form will be sorted.
% Output:
% X (N,n)-matrix containing the ordered Schur-vectors
% columnwise
% RR (N,N) ordered Schur matrix
%
% Written by Bernhard Reuter, Theoretical Physics II,
% University of Kassel, 2017
% -----------------------------------------------------------------------
class_t1 = class(P) ;
class_t2 = class(sd) ;
if (strcmp(class_t1,class_t2))
class_t = class_t1 ;
else
error('do_schur:DataTypeError', ...
['class(P) is not equal class(sd)! This will lead to ', ...
'numeric precision errors!'])
end
function r = num_t(expression)
if (nargin > 0)
if(strcmpi(class_t,'mp')), r = mp(expression) ;
else
if isnumeric(expression)
r = expression ;
else
r = eval(expression) ;
end ;
end ;
else
r = class_t ;
end ;
end % num_t
% -----------------------------------------------------------------------
% assertions
N = size(P,1) ;
assert( size(P,1)==size(P,2), 'do_schur:MatrixShapeError1', ...
'P matrix isnt quadratic!' )
assert( size(sd,1)==size(P,1), 'do_schur:MatrixShapeError2', ...
'sd vector length doesnt match with the shape of P!' )
assert(all(sum(P,2) > num_t('0.0')),'do_schur:ZeroError1', ...
'Not all rows of P are >0!')
assert(all(sd > num_t('0.0')),'do_schur:ZeroError2', ...
'Not all elements of sd are >0!')
% weight the stochastic matrix P_stoch by the initial distribution sd
Pd=diag(sqrt(sd))*P*diag(num_t('1.0')./sqrt(sd)) ;
assert(isa(Pd,num_t),'do_schur:Pd_DataTypeError', ...
'Variable is type %s not %s', class(Pd),num_t)
name=strcat(fileid,'-','Pd','.txt') ;
save_t(name,Pd,'-ascii')
% make a Schur decomposition of Pd
[Q, R]=schur(Pd) ;
assert(isa(Q,num_t),'do_schur:Q_DataTypeError', ...
'Variable is type %s not %s',class(Q),num_t)
assert(isa(R,num_t),'do_schur:R_DataTypeError', ...
'Variable is type %s not %s',class(R),num_t)
name=strcat(fileid,'-','Q','.txt') ;
save_t(name,Q,'-ascii')
name=strcat(fileid,'-','R','.txt') ;
save_t(name,R,'-ascii')
% sort the Schur matrix and vectors
target=num_t('1.0') ;
[ QQ, RR, ap ] = SRSchur_num_t( Q, R, target, b ) ;
assert(isa(QQ,num_t),'do_schur:QQ_DataTypeError', ...
'Variable is type %s not %s',class(QQ),num_t)
assert(isa(RR,num_t),'do_schur:RR_DataTypeError', ...
'Variable is type %s not %s',class(RR),num_t)
assert((any(ap > num_t('1.0')) == 0),'do_schur:ap_ValueError', ...
'ap from SRSchur exceeds one')
name=strcat(fileid,'-','QQ','.txt') ;
save_t(name,QQ,'-ascii')
name=strcat(fileid,'-','RR','.txt') ;
save_t(name,RR,'-ascii')
% if b isn't =0, the Schurmatrix and Schurvectors are only partially
% sorted, therefore one doesn't need the whole Schurvector matrix
if b ~= 0
% mind that b can be <0!
absB = abs(b) ;
% take only the Schurvectors belonging to the sorted part of the
% Schurmatrix
QQ = QQ(:,1:absB) ;
end
% orthonormalize the sorted Schur vectors QQ
% by modified Gram-Schmidt-Orthonormalization
QQ_orthonorm=gram_schmidt_mod(QQ,sd) ;
assert(isa(QQ_orthonorm,num_t),'do_schur:QQ_orthonorm_DataTypeError', ...
'Variable is type %s not %s', class(QQ_orthonorm),num_t)
dummy = ((QQ_orthonorm'*QQ_orthonorm - num_t(eye(size(QQ_orthonorm,2)))) ...
> (num_t('10000')*eps(num_t))) ;
assert(~any(dummy(:)),'do_schur:QQ_orthonorm_OrthogonalityError', ...
'Schurvectors appear to not be orthogonal')
clearvars dummy
name=strcat(fileid,'-','QQ_orthonorm','.txt') ;
save_t(name,QQ_orthonorm,'-ascii')
% transform the orthonormalized Schur vectors of Pd back
% -> orthonormalized Schur vectors X of P_stoch
X=diag(num_t('1.0')./sqrt(sd))*QQ_orthonorm ;
assert( size(X,1)==N, 'do_schur:MatrixShapeError4', ...
'The shape (%d,%d) of X doesnt match with the shape (%d,%d) of P!', ...
size(X,1), size(X,2), N, N)
assert(isa(X,num_t),'do_schur:X_DataTypeError', ...
'X is type %s not %s',class(X),num_t)
dummy = (X'*diag(sd)*X - num_t(eye(size(X,2))) > (num_t('10000')*eps(num_t))) ;
assert(~any(dummy(:)),'do_schur:X_OrthogonalityError', ...
'Schurvectors appear to not be orthogonal')
clearvars dummy
dummy = ( abs(X(:,1) - 1) < ( num_t('1000') * eps(num_t) ) ) ;
assert(all(dummy(:)), 'do_schur:FirstColumnError', ...
'X(:,1) isnt equal 1!')
name=strcat(fileid,'-','X','.txt') ;
save_t(name,X,'-ascii')
end