0001 function checkmanifold(M)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027 assert(isstruct(M), 'M must be a structure.');
0028
0029
0030 list_of_functions = {'name', 'dim', 'inner', 'norm', 'typicaldist', ...
0031 'proj', 'tangent', 'egrad2rgrad', 'retr', ...
0032 'rand', 'randvec', 'zerovec', 'lincomb'};
0033 for k = 1 : numel(list_of_functions)
0034 field = list_of_functions{k};
0035 if ~(isfield(M, field) && isa(M.(field), 'function_handle'))
0036 fprintf('M.%s must be a function handle.\n', field);
0037 end
0038 end
0039
0040
0041 list_of_functions = {'dist', 'ehess2rhess', 'exp', 'log', 'hash', ...
0042 'transp', 'pairmean', 'vec', 'mat', ...
0043 'vecmatareisometries'};
0044 for k = 1 : numel(list_of_functions)
0045 field = list_of_functions{k};
0046 if ~(isfield(M, field) && isa(M.(field), 'function_handle'))
0047 fprintf(['M.%s should ideally (but does not have to) be ' ...
0048 'a function handle.\n'], field);
0049 end
0050 end
0051
0052
0053 try
0054 x = M.rand();
0055 v = M.randvec(x);
0056 fprintf('Random tangent vector norm: %g (should be 1).\n', ...
0057 M.norm(x, v));
0058 z = M.lincomb(x, 1, v, -1, v);
0059 fprintf('norm(v - v)_x = %g (should be 0).\n', M.norm(x, z));
0060 catch up
0061 fprintf('Couldn''t check rand, randvec, lincomb.\n');
0062 end
0063
0064
0065 try
0066 x = M.rand();
0067
0068
0069 u = M.randvec(x);
0070 v = M.randvec(x);
0071 uv = M.inner(x, u, v);
0072 vu = M.inner(x, v, u);
0073 fprintf('<u, v>_x = %g, <v, u>_x = %g, difference = %g (should be 0).\n', uv, vu, uv-vu);
0074
0075
0076 a = randn();
0077 b = randn();
0078 w = M.lincomb(x, a, u, b, v);
0079 z = M.randvec(x);
0080 wz = M.inner(x, w, z);
0081 wzbis = a*M.inner(x, u, z) + b*M.inner(x, v, z);
0082 fprintf('<au+bv, z>_x = %g, a<u, z>_x + b<v, z>_x = %g, difference = %g (should be 0).\n', wz, wzbis, wz-wzbis);
0083
0084
0085
0086
0087 catch up
0088 fprintf('Couldn''t check inner.\n');
0089 end
0090
0091
0092 try
0093 x = M.rand();
0094 v = M.randvec(x);
0095 va = M.tangent2ambient(x, v);
0096 vp = M.proj(x, va);
0097 v_min_vp = M.lincomb(x, 1, v, -1, vp);
0098 df = M.norm(x, v_min_vp);
0099 fprintf('Norm of tangent vector minus its projection to tangent space: %g (should be zero).\n', df);
0100
0101
0102
0103
0104
0105 catch up
0106 fprintf('Couldn''t check tangent2ambient, proj, norm\n');
0107 end
0108
0109
0110 try
0111 x = M.rand();
0112 v = M.randvec(x);
0113 for t = logspace(-8, 1, 10)
0114 y = M.exp(x, v, t);
0115 d = M.dist(x, y);
0116 err = d - abs(t)*M.norm(x, v);
0117 fprintf(['dist(x, M.exp(x, v, t)) - abs(t)*M.norm(x, v) = ' ...
0118 '%g (t = %.1e; should be zero for small enough t).\n'], ...
0119 err, t);
0120 end
0121 catch up
0122 fprintf('Couldn''t check exp and dist.\n');
0123
0124
0125
0126 end
0127
0128
0129 try
0130 x = M.rand();
0131 u = M.randvec(x);
0132 v = M.randvec(x);
0133 U = M.vec(x, u);
0134 V = M.vec(x, v);
0135 if ~iscolumn(U) || ~iscolumn(V)
0136 fprintf('M.vec should return column vectors: they are not.\n');
0137 end
0138 if ~isreal(U) || ~isreal(V)
0139 fprintf('M.vec should return real vectors: they are not real.\n');
0140 end
0141 fprintf(['Unless otherwise stated, M.vec seems to return real ' ...
0142 'column vectors, as intended.\n']);
0143 ru = M.norm(x, M.lincomb(x, 1, M.mat(x, U), -1, u));
0144 rv = M.norm(x, M.lincomb(x, 1, M.mat(x, V), -1, v));
0145 fprintf(['Checking mat/vec are inverse pairs: ' ...
0146 '%g, %g (should be two zeros).\n'], ru, rv);
0147 a = randn(1);
0148 b = randn(1);
0149 fprintf('Checking if vec is linear: %g (should be zero).\n', ...
0150 norm(M.vec(x, M.lincomb(x, a, u, b, v)) - (a*U + b*V)));
0151 if M.vecmatareisometries()
0152 fprintf('M.vecmatareisometries says true.\n');
0153 else
0154 fprintf('M.vecmatareisometries says false.\n');
0155 end
0156 fprintf('If true, this should be zero: %g.\n', ...
0157 U(:).'*V(:) - M.inner(x, u, v));
0158 catch up
0159 fprintf('Couldn''t check mat, vec, vecmatareisometries.\n');
0160 end
0161
0162
0163 dim_threshold = 200;
0164 if M.dim() <= dim_threshold
0165 x = M.rand();
0166 n = M.dim() + 1;
0167 B = cell(n, 1);
0168 for k = 1 : n
0169 B{k} = M.randvec(x);
0170 end
0171 G = grammatrix(M, x, B);
0172 eigG = sort(real(eig(G)), 'descend');
0173 fprintf('Testing M.dim() (works best when dimension is small):\n');
0174 fprintf('\tIf this number is machine-precision zero, then M.dim() may be too large: %g\n', eigG(n-1));
0175 fprintf('\tIf this number is not machine-precision zero, then M.dim() may be too small: %g\n', eigG(n));
0176 else
0177 fprintf('M.dim() not tested because it is > %d.\n', dim_threshold);
0178 end
0179
0180
0181 fprintf('It is recommended also to call checkretraction.\n');
0182
0183 end
0184