1function test_spqr_coverage 2%TEST_SPQR_COVERAGE statement coverage test of spqr_rank functions 3% Example 4% test_spqr_coverage 5% 6% This test exercises all but a handful of statements in the spqr_rank toolbox. 7% A lot of output is generated, and some of it will contain error messages. 8% This is expected, since the spqr_* functions do a lot of error checking, and 9% this code exercises essentially all of those tests with many pathological 10% matrices. Sample output of this test is in private/test_spqr_coverage.txt. 11% 12% The only statements not tested are a few marked "untested" in the code. 13% These are for a few rare error conditions which might occur, or which might 14% actually be theoretically impossible. Since we have no proof that these 15% conditions cannot occur, we include code to check for them and handle the 16% potential errors. But we cannot test this error-handling code itself. 17% 18% After running this test, you can view the Coverage reports by going to the 19% Current Folder window, selecting the Gear symbol and selecting Reports -> 20% Coverage Report. Some files report less than 100%, but this is the fault of 21% MATLAB not accounting for unreachable statements. For example, code such as 22% this reports the "end" statement as untested: 23% 24% if (...) 25% return 26% end 27% 28% Coverage is close to 100% for all files in the toolbox except for the testing/ 29% demonstration files (test_spqr_coverage, demo_spqr_rank and test_spqr_rank). 30% 31% See also demo_spqr_rank, test_spqr_rank 32 33% Copyright 2012, Leslie Foster and Timothy A Davis. 34 35profile clear 36profile on 37 38fprintf (['Exercises all but a handful of statements in the ', ... 39 'spqr_rank toolbox\n']) ; 40 41nfail = 0 ; 42 43opts = struct('figures', 0, ... % no figures 44 'null_spaces', 1, ... 45 'doprint', -1, ... % absolutely no printing 46 'start_with_A_transpose', 0, ... 47 'implicit_null_space_basis', 1, ... 48 'repeatable', 1 , ... 49 'nsvals', 1 , ... 50 'get_details', 0, ... 51 'tol_norm_type', 2) ; 52nfail = nfail + demo_spqr_rank (3, opts) ; 53opts.tol_norm_type = 1 ; 54spqr_rank_opts (opts) ; % exercise spqr_rank_opts with tol_norm_type = 1 55nfail = nfail + demo_spqr_rank (3, opts) ; 56 57A = magic (5) ; 58A = A (1:4,:) ; 59b = (1:4)' ; 60 61try 62 % incorrect usage: "A" not defined 63 x = spqr_basic ; %#ok 64 fprintf ('usage error not caught\n') ; 65 nfail = nfail + 1 ; 66catch me 67 fprintf ('Expected error caught: %s\n', me.message) ; 68end 69 70try 71 % broken stats struct 72 stats.flag = 999 ; 73 spqr_rank_stats (stats,1) ; 74 nfail = nfail + 1 ; 75catch me 76 fprintf ('Expected error caught: %s\n', me.message) ; 77end 78 79try 80 % incorrect usage: too many arguments 81 x = spqr_basic (A, b, 1, opts) ; %#ok 82 fprintf ('usage error not caught\n') ; 83 nfail = nfail + 1 ; 84catch me 85 fprintf ('Expected error caught: %s\n', me.message) ; 86end 87 88try 89 % incorrect usage: A and b mismatch 90 c = (1:5)' ; 91 x = spqr_basic (A, c) ; %#ok 92 fprintf ('usage error not caught\n') ; 93 nfail = nfail + 1 ; 94catch me 95 fprintf ('Expected error caught: %s\n', me.message) ; 96end 97 98try 99 % incorrect usage: 3rd argument not scalar 100 x = spqr_basic (A, b, [1 2]) ; %#ok 101 fprintf ('usage error not caught\n') ; 102 nfail = nfail + 1 ; 103catch me 104 fprintf ('Expected error caught: %s\n', me.message) ; 105end 106 107try 108 % incorrect usage: too many arguments 109 [U, S, V, stats] = spqr_ssi (A, 1, opts) ; %#ok 110 fprintf ('usage error not caught\n') ; 111 nfail = nfail + 1 ; 112catch me 113 fprintf ('Expected error caught: %s\n', me.message) ; 114end 115 116try 117 % incorrect usage: too many arguments 118 x = spqr_cod (A, b, 1, opts) ; %#ok 119 fprintf ('usage error not caught\n') ; 120 nfail = nfail + 1 ; 121catch me 122 fprintf ('Expected error caught: %s\n', me.message) ; 123end 124 125try 126 % incorrect usage: too many arguments 127 N = spqr_null (A, b, opts) ; %#ok 128 fprintf ('usage error not caught\n') ; 129 nfail = nfail + 1 ; 130catch me 131 fprintf ('Expected error caught: %s\n', me.message) ; 132end 133 134try 135 % incorrect usage: too many arguments 136 [U,S,V,stats] = spqr_ssp (A, [ ], 4, pi, opts) ; %#ok 137 fprintf ('usage error not caught\n') ; 138 nfail = nfail + 1 ; 139catch me 140 fprintf ('Expected error caught: %s\n', me.message) ; 141end 142 143try 144 % incorrect usage: too many arguments 145 x = spqr_pinv (A, b, 1, opts) ; %#ok 146 fprintf ('usage error not caught\n') ; 147 nfail = nfail + 1 ; 148catch me 149 fprintf ('Expected error caught: %s\n', me.message) ; 150end 151 152try 153 % incorrect usage: A must be square 154 [U, S, V, stats] = spqr_ssi (A) ; %#ok 155 fprintf ('usage error not caught\n') ; 156 nfail = nfail + 1 ; 157catch me 158 fprintf ('Expected error caught: %s\n', me.message) ; 159end 160 161%------------------------------------------------------------------------------- 162% check basic correct usage 163%------------------------------------------------------------------------------- 164 165x = spqr_basic (A, b, [ ]) ; 166if (norm (A*x-b) > 1e-13) 167 nfail = nfail + 1 ; 168end 169 170x = spqr_basic (A, b, eps) ; 171if (norm (A*x-b) > 1e-13) 172 nfail = nfail + 1 ; 173end 174 175x = spqr_basic (A, b, struct ('tol_norm_type', 1)) ; 176if (norm (A*x-b) > 1e-13) 177 nfail = nfail + 1 ; 178end 179 180c = spqr_ssp (A,1) ; 181if (abs (c - norm (A)) > 0.01) 182 nfail = nfail + 1 ; 183end 184 185opts2.get_details = 1 ; 186fprintf ('\nopts for spqr_ssi with get_details = 1:\n') ; 187spqr_rank_opts (opts2) ; 188R = qr (A (:, 1:4)) ; 189[U, S, V, stats] = spqr_ssi (R, opts2) ; %#ok 190c = stats.normest_R ; 191spqr_rank_stats (stats) ; 192if (abs (c - norm (R)) / norm (R) > 0.01) 193 nfail = nfail + 1 ; 194end 195 196[s, stats] = spqr_ssi (inf) ; %#ok 197 198for t = 0:2 199 opts2.get_details = t ; 200 opts2.tol = 1e-6 ; 201 fprintf ('\nopts for spqr_cod with get_details = %d:\n', t) ; 202 spqr_rank_opts (opts2) ; 203 [x, stats] = spqr_cod (A, b, opts2) ; 204 e = norm (A*x-b) ; 205 fprintf ('\nresults from spqr_cod with get_details = %d: err %g\n', t, e) ; 206 spqr_rank_stats (stats) ; 207 if (e > 1e-12) 208 nfail = nfail + 1 ; 209 end 210 211 [x stats N NT] = spqr_pinv (A, b, opts2) ; %#ok 212 e = norm (A*x-b) ; 213 fprintf ('\nresults from spqr_pinv with get_details = %d: err %g\n', t, e) ; 214 spqr_rank_stats (stats, 1) ; 215 if (e > 1e-12) 216 nfail = nfail + 1 ; 217 end 218 219 [x stats NT] = spqr_basic (A, b, opts2) ; %#ok 220 e = norm (A*x-b) ; 221 fprintf ('\nresults from spqr_basic with get_details = %d: err %g\n', t, e); 222 spqr_rank_stats (stats, 1) ; 223 if (e > 1e-12) 224 nfail = nfail + 1 ; 225 end 226 227end 228 229C = magic (4) ; % matrix of rank 3 230opts3.implicit_null_space_basis = 0 ; 231for t = 0:2 232 opts3.get_details = t ; 233 [N stats] = spqr_null (C, opts3) ; %#ok 234 if (norm (C*N) > 1e-12) 235 nfail = nfail + 1 ; 236 end 237end 238 239[s stats] = spqr_ssp (A, [ ], 0) ; %#ok 240s = spqr_ssp (A, [ ], 0) ; %#ok 241s = spqr_ssp (sparse (pi)) ; 242if (abs (s - pi) > 1e-12) 243 nfail = nfail + 1 ; 244end 245 246s = spqr_ssp (A, 10) ; %#ok 247 248C = magic (4) ; % matrix of rank 3 249B = (1:4) ; 250opts6.implicit_null_space_basis = 0 ; 251opts6.start_with_A_transpose = 1 ; 252opts6.ssi_tol = 1e-6 ; 253opts6.repeatable = 0 ; 254opts6.get_details = 1 ; 255fprintf ('\nopts for spqr_null with explicit basis N:\n') ; 256spqr_rank_opts (opts6) ; 257N1 = spqr_null (C) ; 258N2 = spqr_null (C, opts6) ; 259e = 0 ; 260e = max (e, norm (N2'*C - spqr_null_mult (N1, C))) ; 261e = max (e, norm (N2'*C - spqr_null_mult (N1, C, 0))) ; 262e = max (e, norm (N2*B - spqr_null_mult (N1, B, 1))) ; 263e = max (e, norm (B'*N2' - spqr_null_mult (N1, B', 2))) ; 264e = max (e, norm (C*N2 - spqr_null_mult (N1, C, 3))) ; 265e = max (e, norm (B*N2 - spqr_null_mult (N1, B, 3))) ; 266 267e = max (e, norm (N2'*C - spqr_null_mult (N2, C))) ; 268e = max (e, norm (N2'*C - spqr_null_mult (N2, C, 0))) ; 269e = max (e, norm (N2*B - spqr_null_mult (N2, B, 1))) ; 270e = max (e, norm (B'*N2' - spqr_null_mult (N2, B', 2))) ; 271e = max (e, norm (C*N2 - spqr_null_mult (N2, C, 3))) ; 272e = max (e, norm (B*N2 - spqr_null_mult (N2, B, 3))) ; 273 274if (e > 1e-12) 275 nfail = nfail + 1 ; 276end 277 278try 279 % unrecognized method 280 x = spqr_null_mult (N1, C, 42) ; %#ok 281catch me 282 fprintf ('Expected error caught: %s\n', me.message) ; 283end 284 285try 286 % unrecognized method 287 x = spqr_null_mult (N2, C, 42) ; %#ok 288catch me 289 fprintf ('Expected error caught: %s\n', me.message) ; 290end 291 292Ngunk = N1 ; 293Ngunk.kind = 'gunk' ; 294 295try 296 % unrecognized struct 297 x = spqr_null_mult (Ngunk, C, 0) ; %#ok 298catch me 299 fprintf ('Expected error caught: %s\n', me.message) ; 300end 301 302try 303 % unrecognized struct 304 x = spqr_null_mult (Ngunk, B, 1) ; %#ok 305catch me 306 fprintf ('Expected error caught: %s\n', me.message) ; 307end 308 309try 310 % unrecognized struct 311 x = spqr_null_mult (Ngunk, B', 2) ; %#ok 312catch me 313 fprintf ('Expected error caught: %s\n', me.message) ; 314end 315 316try 317 % unrecognized struct 318 x = spqr_null_mult (Ngunk, C, 3) ; %#ok 319catch me 320 fprintf ('Expected error caught: %s\n', me.message) ; 321end 322 323C = magic (4) ; 324C (:,1) = 0.2 * C (:,2) + pi * C (:,3) ; 325b = (1:4)' ; 326opts6.get_details = 2 ; 327[x,stats,N1,NT1] = spqr_cod (C,b) ; %#ok 328[x,stats,N2,NT2] = spqr_cod (C,b, opts6) ; %#ok 329fprintf ('\ndetailed stats from spqr_cod:\n') ; 330spqr_rank_stats (stats, 1) ; 331e = 0 ; 332 333e = max (e, norm (N2'*C - spqr_null_mult (N1, C))) ; 334e = max (e, norm (N2'*C - spqr_null_mult (N1, C, 0))) ; 335e = max (e, norm (N2*B - spqr_null_mult (N1, B, 1))) ; 336e = max (e, norm (B'*N2' - spqr_null_mult (N1, B', 2))) ; 337e = max (e, norm (C*N2 - spqr_null_mult (N1, C, 3))) ; 338e = max (e, norm (B*N2 - spqr_null_mult (N1, B, 3))) ; 339 340e = max (e, norm (N2'*C - spqr_null_mult (N2, C))) ; 341e = max (e, norm (N2'*C - spqr_null_mult (N2, C, 0))) ; 342e = max (e, norm (N2*B - spqr_null_mult (N2, B, 1))) ; 343e = max (e, norm (B'*N2' - spqr_null_mult (N2, B', 2))) ; 344e = max (e, norm (C*N2 - spqr_null_mult (N2, C, 3))) ; 345e = max (e, norm (B*N2 - spqr_null_mult (N2, B, 3))) ; 346 347e = max (e, norm (NT2'*C - spqr_null_mult (NT1, C))) ; 348e = max (e, norm (NT2'*C - spqr_null_mult (NT1, C, 0))) ; 349e = max (e, norm (NT2*B - spqr_null_mult (NT1, B, 1))) ; 350e = max (e, norm (B'*NT2' - spqr_null_mult (NT1, B', 2))) ; 351e = max (e, norm (C*NT2 - spqr_null_mult (NT1, C, 3))) ; 352e = max (e, norm (B*NT2 - spqr_null_mult (NT1, B, 3))) ; 353 354Nsparse = spqr_explicit_basis (N1) ; 355NTsparse = spqr_explicit_basis (NT1) ; 356 357e = max (e, norm (N2'*C - spqr_null_mult (Nsparse, C))) ; 358e = max (e, norm (N2'*C - spqr_null_mult (Nsparse, C, 0))) ; 359e = max (e, norm (N2*B - spqr_null_mult (Nsparse, B, 1))) ; 360e = max (e, norm (B'*N2' - spqr_null_mult (Nsparse, B', 2))) ; 361e = max (e, norm (C*N2 - spqr_null_mult (Nsparse, C, 3))) ; 362e = max (e, norm (B*N2 - spqr_null_mult (Nsparse, B, 3))) ; 363 364e = max (e, norm (NT2'*C - spqr_null_mult (NTsparse, C))) ; 365e = max (e, norm (NT2'*C - spqr_null_mult (NTsparse, C, 0))) ; 366e = max (e, norm (NT2*B - spqr_null_mult (NTsparse, B, 1))) ; 367e = max (e, norm (B'*NT2' - spqr_null_mult (NTsparse, B', 2))) ; 368e = max (e, norm (C*NT2 - spqr_null_mult (NTsparse, C, 3))) ; 369e = max (e, norm (B*NT2 - spqr_null_mult (NTsparse, B, 3))) ; 370 371Nfull = spqr_explicit_basis (N1, 'full') ; 372NTfull = spqr_explicit_basis (NT1, 'full') ; 373 374e = max (e, norm (N2'*C - spqr_null_mult (Nfull, C))) ; 375e = max (e, norm (N2'*C - spqr_null_mult (Nfull, C, 0))) ; 376e = max (e, norm (N2*B - spqr_null_mult (Nfull, B, 1))) ; 377e = max (e, norm (B'*N2' - spqr_null_mult (Nfull, B', 2))) ; 378e = max (e, norm (C*N2 - spqr_null_mult (Nfull, C, 3))) ; 379e = max (e, norm (B*N2 - spqr_null_mult (Nfull, B, 3))) ; 380 381e = max (e, norm (NT2'*C - spqr_null_mult (NTfull, C))) ; 382e = max (e, norm (NT2'*C - spqr_null_mult (NTfull, C, 0))) ; 383e = max (e, norm (NT2*B - spqr_null_mult (NTfull, B, 1))) ; 384e = max (e, norm (B'*NT2' - spqr_null_mult (NTfull, B', 2))) ; 385e = max (e, norm (C*NT2 - spqr_null_mult (NTfull, C, 3))) ; 386e = max (e, norm (B*NT2 - spqr_null_mult (NTfull, B, 3))) ; 387 388Nfull = spqr_explicit_basis (Nfull, 'full') ; % doesn't change Nfull 389e = max (e, norm (N2'*C - spqr_null_mult (Nfull, C))) ; 390 391if (e > 1e-12) 392 nfail = nfail + 1 ; 393end 394 395% test spqr_wrapper 396A = sparse (magic (4)) ; 397[Q,R,C,p,info] = spqr_wrapper (A, [ ], eps, 'discard', 0) ; %#ok 398if (~isequal (size (C), [4 0]) || norm (R'*R - A'*A, 1) > 1e-12) 399 nfail = nfail + 1 ; 400end 401 402fprintf ('\ndefault ') ; 403spqr_rank_opts 404opts = spqr_rank_opts (struct, 1) ; 405fprintf ('\n(default for spqr_ssp): ') ; 406spqr_rank_opts (opts) ; 407 408fprintf ('\ndescription of statistics:\n') ; 409spqr_rank_stats 410spqr_rank_stats ('ssi') ; 411spqr_rank_stats ('ssp') ; 412spqr_rank_stats ('null') ; 413spqr_rank_stats (1) ; 414spqr_rank_stats (2) ; 415 416fprintf ('\n---------------------------------------------------------------\n'); 417fprintf ('test that illustrates the rare case of a miscalculated rank:\n') ; 418install_SJget ; 419Problem = SJget (182) ; 420display (Problem) 421A = Problem.A ; 422opts8.get_details = 1; 423opts8.tol = max(size(A))*eps(Problem.svals(1)); 424[N,stats] = spqr_null (A,opts8) ; %#ok 425opts8.repeatable = 0 ; 426[N,stats] = spqr_null (A,opts8) ; %#ok 427opts8.repeatable = 384 ; 428[N,stats] = spqr_null (A,opts8) ; %#ok 429spqr_rank_stats (stats,1) ; 430if (stats.flag == 1) 431 rank_svd = SJrank (Problem,stats.tol_alt) ; 432else 433 rank_svd = SJrank (Problem,stats.tol) ; 434end 435if stats.rank ~= rank_svd 436 fprintf ('\nexpected rank mismatch %d %d\n', stats.rank, rank_svd) ; 437end 438 439% another rare case 440fprintf ('\n---------------------------------------------------------------\n'); 441fprintf ('another rare case:\n') ; 442Problem = SJget (240) ; 443display (Problem) ; 444A = Problem.A ; 445m = size (A,1) ; 446b = ones (m, 1) ; 447opts9.repeatable = 22 ; 448opts9.get_details = 1 ; 449[x stats] = spqr_cod (A, b, opts9) ; %#ok 450spqr_rank_stats (stats,1) ; 451 452% another rare case: early return from spqr_ssi 453fprintf ('\n---------------------------------------------------------------\n'); 454fprintf ('another rare case:\n') ; 455Problem = SJget (242) ; 456display (Problem) ; 457A = Problem.A ; 458m = size (A,1) ; 459b = ones (m, 1) ; 460opts9.repeatable = 1 ; 461[N stats] = spqr_null (A, opts9) ; %#ok 462spqr_rank_stats (stats,1) ; 463x = spqr_pinv (A', b, opts9) ; %#ok 464[Q,R,C,p,info] = spqr_wrapper (A', [ ], eps, 'discard Q', 2) ; %#ok 465rank_spqr = size (R,1) ; 466R11 = R (:, 1:rank_spqr) ; 467opts9.tol = 1e-14 ; 468[s stats] = spqr_ssi (R11, opts9) ; %#ok 469spqr_rank_stats (stats,1) ; 470 471% test rare case in spqr_ssi 472fprintf ('\n---------------------------------------------------------------\n'); 473fprintf ('Near-overflow in spqr_ssi, 2nd test:\n') ; 474Problem = SJget (137) ; 475display (Problem) ; 476opts10.repeatable = 1 ; 477opts10.get_details = 1 ; 478A = Problem.A ; 479[m n] = size (A) ; 480tol = max(m,n) * eps (normest (A,0.01)) ; 481for t = 0:1 482 [Q,R,C,p,info] = spqr_wrapper (A, [ ], tol, 'discard Q', 2) ; %#ok 483 k = size (R,1) ; 484 R = R (:, 1:k) ; 485 [s, stats] = spqr_ssi (R, opts10) ; %#ok 486 spqr_rank_stats (stats,1) ; 487 [s, stats] = spqr_ssi (R, opts10) ; %#ok 488 spqr_rank_stats (stats,1) ; 489 A = A' ; 490end 491 492% test flag=1 case: rank correct with an alternate tolerance 493fprintf ('\n---------------------------------------------------------------\n'); 494fprintf ('case where rank appears correct but with alternate tolerance:\n') ; 495Problem = SJget (183) ; 496display (Problem) ; 497A = Problem.A ; 498b = ones (size (A,1), 1) ; 499opts11.get_details = 1 ; 500opts11.repeatable = 3 ; 501[x stats] = spqr_cod (A,b, opts11) ; %#ok 502spqr_rank_stats (stats,1) ; 503 504if (nfail > 0) 505 error ('nfail: %d\n', nfail) ; 506end 507 508%------------------------------------------------------------------------------- 509% exhaustive tests on a handful of matrices 510%------------------------------------------------------------------------------- 511 512index = SJget ; 513 514% these matrices trigger lots of special cases in the spqr_rank toolbox: 515% 229 Regtools/foxgood_100 516% 241 Regtools/heat_100 517% 245 Regtools/i_laplace_100 518% 249 Regtools/parallax_100 519% 242 Regtools/heat_200 520% 173 Sandia/oscil_dcop_24 521% 183 Sandia/oscil_dcop_34 522% 263 Regtools/wing_500 523idlist = [ 215 229 241 245 249 242 173 183 263 134 ] ; 524 525fprintf ('\nPlease wait ...\n') ; 526for k = 1:length (idlist) 527 id = idlist (k) ; 528 fprintf ('%2d of %2d : id: %4d %s/%s\n', k, length (idlist), ... 529 id, index.Group {id}, index.Name {id}) ; 530 nfail = nfail + test_spqr_rank (id, 0) ; 531end 532 533if (nfail > 0) 534 error ('One or more tests failed!') ; 535end 536 537fprintf ('\nAll tests passed.\n') ; 538 539profile off 540 541%------------------------------------------------------------------------------- 542% coverage results 543%------------------------------------------------------------------------------- 544 545% This test_spqr_coverage script obtains near-100% coverage. 546% 547% These files obtain 100% coverage: 548% 549% spqr_rank_opts.m 550% spqr_explicit_basis.m 551% private/tol_is_default 552% private/spqr_rank_order_fields 553% private/spqr_rank_form_basis 554% private/spqr_rank_deflation 555% private/spqr_rank_assign_stats 556% private/spqr_wrapper 557 558% These files are effectively 100%. The only non-tested statements are "end" 559% statements that appear immediate after 'error' or 'break' statements, which 560% are not reachable. MATLAB should report 100% coverage in this case, but it 561% doesn't. This is a failure of MATLAB, not our code... 562% 563% spqr_ssp 564% spqr_null 565% spqr_null_mult 566% spqr_basic 567% spqr_rank_stats 568% private/spqr_rank_get_inputs 569 570% These files are not yet 100%, and "cannot" be. These codes have a few checks 571% for errors that should "never" occur. We know of no matrices that trigger 572% these conditions, but also no proof that the conditions cannot occur. Lines 573% that are untested are marked with "% untested" comments. 574% 575% spqr_cod 576% spqr_ssi 577% spqr_pinv 578 579% One part of this file is not tested (error handling if 'savepath' fails): 580% 581% private/install_SJget 582 583% These files are not fully covered, but do not need to be (test/demo code): 584% 585% test_spqr_coverage 586% test_spqr_rank 587% demo_spqr_rank 588% files in the SJget folder 589 590