Mercurial > octave

annotate scripts/sparse/sprandsym.m @ 33625:d213a148b3f1 default tip @

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ensure exp. terminal widget has focus at startup * main-window.cc (main_window): call focus_command_window only if event loop is idle by using a single shot timer * main-window.h: make focus_command_window a public slot
author Torsten Lilge <ttl-octave@mailbox.org>
date Sun, 26 May 2024 02:29:44 +0200
parents 2e484f9f1f18
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
1 ########################################################################
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
2 ##
32632
2e484f9f1f18 maint: update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 31706
diff changeset
3 ## Copyright (C) 2004-2024 The Octave Project Developers
27918
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
4 ##
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
5 ## See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
6 ## distribution or <https://octave.org/copyright/>.
5610
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman]
dbateman
parents:
diff changeset
7 ##
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6498
diff changeset
8 ## This file is part of Octave.
5610
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman]
dbateman
parents:
diff changeset
9 ##
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
10 ## Octave is free software: you can redistribute it and/or modify it
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6498
diff changeset
11 ## under the terms of the GNU General Public License as published by
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
12 ## the Free Software Foundation, either version 3 of the License, or
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
13 ## (at your option) any later version.
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6498
diff changeset
14 ##
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6498
diff changeset
15 ## Octave is distributed in the hope that it will be useful, but
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6498
diff changeset
16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22323
diff changeset
18 ## GNU General Public License for more details.
5610
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman]
dbateman
parents:
diff changeset
19 ##
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman]
dbateman
parents:
diff changeset
20 ## You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6498
diff changeset
21 ## along with Octave; see the file COPYING. If not, see
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23220
diff changeset
22 ## <https://www.gnu.org/licenses/>.
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
23 ##
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
24 ########################################################################
5610
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman]
dbateman
parents:
diff changeset
25
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman]
dbateman
parents:
diff changeset
26 ## -*- texinfo -*-
30875
5d3faba0342e doc: Ensure documentation lists output argument when it exists for all m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
27 ## @deftypefn {} {@var{S} =} sprandsym (@var{n}, @var{d})
5d3faba0342e doc: Ensure documentation lists output argument when it exists for all m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
28 ## @deftypefnx {} {@var{S} =} sprandsym (@var{s})
18809
53af80da6781 doc: Update documentation of sparse functions including seealso links.
Rik <rik@octave.org>
parents: 17744
diff changeset
29 ## Generate a symmetric random sparse matrix.
5610
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman]
dbateman
parents:
diff changeset
30 ##
18809
53af80da6781 doc: Update documentation of sparse functions including seealso links.
Rik <rik@octave.org>
parents: 17744
diff changeset
31 ## The size of the matrix will be @var{n}x@var{n}, with a density of values
53af80da6781 doc: Update documentation of sparse functions including seealso links.
Rik <rik@octave.org>
parents: 17744
diff changeset
32 ## given by @var{d}. @var{d} must be between 0 and 1 inclusive. Values will
53af80da6781 doc: Update documentation of sparse functions including seealso links.
Rik <rik@octave.org>
parents: 17744
diff changeset
33 ## be normally distributed with a mean of zero and a variance of 1.
53af80da6781 doc: Update documentation of sparse functions including seealso links.
Rik <rik@octave.org>
parents: 17744
diff changeset
34 ##
53af80da6781 doc: Update documentation of sparse functions including seealso links.
Rik <rik@octave.org>
parents: 17744
diff changeset
35 ## If called with a single matrix argument, a random sparse matrix is generated
20164
df437a52bcaf doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents: 19833
diff changeset
36 ## wherever the matrix @var{s} is nonzero in its lower triangular part.
18809
53af80da6781 doc: Update documentation of sparse functions including seealso links.
Rik <rik@octave.org>
parents: 17744
diff changeset
37 ## @seealso{sprand, sprandn, spones, sparse}
5610
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman]
dbateman
parents:
diff changeset
38 ## @end deftypefn
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman]
dbateman
parents:
diff changeset
39
6252
738c97e101eb [project @ 2007-01-24 19:58:46 by jwe]
jwe
parents: 6046
diff changeset
40 function S = sprandsym (n, d)
8506
bc982528de11 comment style fixes
John W. Eaton <jwe@octave.org>
parents: 7505
diff changeset
41
28891
de5f2f9a64ff maint: Use same coding style when checking for a minimum of 1 input.
Rik <rik@octave.org>
parents: 28059
diff changeset
42 if (nargin < 1)
6046
34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe]
jwe
parents: 5642
diff changeset
43 print_usage ();
5610
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman]
dbateman
parents:
diff changeset
44 endif
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
45
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
46 if (nargin == 1)
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
47 [i, j] = find (tril (n));
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
48 [nr, nc] = size (n);
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
49 S = sparse (i, j, randn (size (i)), nr, nc);
20231
83792dd9bcc1 Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents: 20164
diff changeset
50 S += tril (S, -1)';
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
51 return;
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
52 endif
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
53
28059
e9d57f6d6353 allow sprand and sprandsym to create empty sparse matrices
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
54 if (!(isscalar (n) && n == fix (n) && n >= 0))
e9d57f6d6353 allow sprand and sprandsym to create empty sparse matrices
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
55 error ("sprandsym: N must be a non-negative integer 0");
e9d57f6d6353 allow sprand and sprandsym to create empty sparse matrices
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
56 endif
e9d57f6d6353 allow sprand and sprandsym to create empty sparse matrices
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
57
e9d57f6d6353 allow sprand and sprandsym to create empty sparse matrices
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
58 if (n == 0)
e9d57f6d6353 allow sprand and sprandsym to create empty sparse matrices
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
59 S = sparse (n, n);
e9d57f6d6353 allow sprand and sprandsym to create empty sparse matrices
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
60 return;
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
61 endif
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
62
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
63 if (d < 0 || d > 1)
13240
32980cbf2338 Use correct function name in error message in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13220
diff changeset
64 error ("sprandsym: density D must be between 0 and 1");
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
65 endif
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
66
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
67 ## Actual number of nonzero entries
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
68 k = round (n^2*d);
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
69
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
70 ## Diagonal nonzero entries, same parity as k
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
71 r = pick_rand_diag (n, k);
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
72
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
73 ## Off diagonal nonzero entries
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
74 m = (k - r)/2;
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
75
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
76 ondiag = randperm (n, r);
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
77 offdiag = randperm (n*(n - 1)/2, m);
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
78
13212
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
79 ## Row index
13213
544304a09e42 Fix offbyones and typos in sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13212
diff changeset
80 i = lookup (cumsum (0:n), offdiag - 1) + 1;
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
81
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
82 ## Column index
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
83 j = offdiag - (i - 1).*(i - 2)/2;
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
84
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
85 diagvals = randn (1, r);
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
86 offdiagvals = randn (1, m);
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
87
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
88 S = sparse ([ondiag, i, j], [ondiag, j, i],
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
89 [diagvals, offdiagvals, offdiagvals], n, n);
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
90
5610
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman]
dbateman
parents:
diff changeset
91 endfunction
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
92
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
93 function r = pick_rand_diag (n, k)
30379
363fb10055df maint: Style check m-files ahead of 7.1 release.
Rik <rik@octave.org>
parents: 29359
diff changeset
94
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
95 ## Pick a random number R of entries for the diagonal of a sparse NxN
13857
9f28f0d05473 sprandsym.m: trivial doc fix
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13240
diff changeset
96 ## symmetric square matrix with exactly K nonzero entries, ensuring
9f28f0d05473 sprandsym.m: trivial doc fix
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13240
diff changeset
97 ## that this R is chosen uniformly over all such matrices.
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
98 ##
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
99 ## Let D be the number of diagonal entries and M the number of
21751
b571fc85953f maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents: 21602
diff changeset
100 ## off-diagonal entries. Then K = D + 2*M. Let A = N*(N-1)/2 be the
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
101 ## number of available entries in the upper triangle of the matrix.
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
102 ## Then, by a simple counting argument, there is a total of
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
103 ##
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
104 ## T = nchoosek (N, D) * nchoosek (A, M)
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
105 ##
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
106 ## symmetric NxN matrices with a total of K nonzero entries and D on
21751
b571fc85953f maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents: 21602
diff changeset
107 ## the diagonal. Letting D range from mod (K,2) through min (N,K), and
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
108 ## dividing by this sum, we obtain the probability P for D to be each
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
109 ## of those values.
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
110 ##
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
111 ## However, we cannot use this form for computation, as the binomial
21751
b571fc85953f maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents: 21602
diff changeset
112 ## coefficients become unmanageably large. Instead, we use the
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
113 ## successive quotients Q(i) = T(i+1)/T(i), which we easily compute to
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
114 ## be
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
115 ##
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
116 ## (N - D)*(N - D - 1)*M
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
117 ## Q = -------------------------------
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
118 ## (D + 2)*(D + 1)*(A - M + 1)
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
119 ##
13212
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
120 ## Then, after prepending 1, the cumprod of these quotients is
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
121 ##
13212
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
122 ## C = [ T(1)/T(1), T(2)/T(1), T(3)/T(1), ..., T(N)/T(1) ]
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
123 ##
13212
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
124 ## Their sum is thus S = sum (T)/T(1), and then C(i)/S is the desired
21751
b571fc85953f maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents: 21602
diff changeset
125 ## probability P(i) for i=1:N. The cumsum will finally give the
13212
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
126 ## distribution function for computing the random number of entries on
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
127 ## the diagonal R.
13220
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
128 ##
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
129 ## Thanks to Zsbán Ambrus <ambrus@math.bme.hu> for most of the ideas
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
130 ## of the implementation here, especially how to do the computation
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
131 ## numerically to avoid overflow.
13212
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
132
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
133 ## Degenerate case
14552
86854d032a37 maint: miscellaneous style fixes for .m files
John W. Eaton <jwe@octave.org>
parents: 14373
diff changeset
134 if (k == 1)
13212
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
135 r = 1;
17312
088d014a7fe2 Use semicolon after "return" statement in core m-files.
Rik <rik@octave.org>
parents: 16037
diff changeset
136 return;
13212
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
137 endif
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
138
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
139 ## Compute the stuff described above
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
140 a = n*(n - 1)/2;
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
141 d = [mod(k,2):2:min(n,k)-2];
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
142 m = (k - d)/2;
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
143 q = (n - d).*(n - d - 1).*m ./ (d + 2)./(d + 1)./(a - m + 1);
13220
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
144
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
145 ## Slight modification from discussion above: pivot around the max in
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
146 ## order to avoid overflow (underflow is fine, just means effectively
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
147 ## zero probabilities).
21602
47c76e56a1de maint: Remove extra space between end of statement and semicolon.
Rik <rik@octave.org>
parents: 20852
diff changeset
148 [~, midx] = max (cumsum (log (q)));
20735
418ae0cb752f Replace ++,-- with in-place operators for performance.
Rik <rik@octave.org>
parents: 20231
diff changeset
149 midx += 1;
13220
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
150 lc = fliplr (cumprod (1./q(midx-1:-1:1)));
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
151 rc = cumprod (q(midx:end));
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
152
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
153 ## Now c = t(i)/t(midx), so c > 1 == [].
410de573089a Avoid overflow in sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13213
diff changeset
154 c = [lc, 1, rc];
13212
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
155 s = sum (c);
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
156 p = c/s;
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
157
13212
e81b93284605 Various minor stylistic improvements to sprandsym.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13208
diff changeset
158 ## Add final d
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
159 d(end+1) = d(end) + 2;
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
160
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
161 ## Pick a random r using this distribution
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
162 r = d(sum (cumsum (p) < rand) + 1);
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
163
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
164 endfunction
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
165
14363
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
166
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
167 %!test
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
168 %! s = sprandsym (10, 0.1);
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
169 %! assert (issparse (s));
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
170 %! assert (issymmetric (s));
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
171 %! assert (size (s), [10, 10]);
13205
abf1e00111dd Completely new implementation of sprandsym
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents: 13064
diff changeset
172 %! assert (nnz (s) / numel (s), 0.1, .01);
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
173
19833
9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
174 ## Test 1-input calling form
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
175 %!test
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
176 %! s = sprandsym (sparse ([1 2 3], [3 2 3], [2 2 2]));
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
177 %! [i, j] = find (s);
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
178 %! assert (sort (i), [2 3]');
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
179 %! assert (sort (j), [2 3]');
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
180
28059
e9d57f6d6353 allow sprand and sprandsym to create empty sparse matrices
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
181 ## Test empty array creation
e9d57f6d6353 allow sprand and sprandsym to create empty sparse matrices
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
182 %!assert (size (sprandsym (0, 0.5)), [0, 0])
e9d57f6d6353 allow sprand and sprandsym to create empty sparse matrices
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
183
19833
9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents: 19697
diff changeset
184 ## Test input validation
28896
90fea9cc9caa test: Add expected error message <Invalid call> to BIST tests for nargin.
Rik <rik@octave.org>
parents: 28891
diff changeset
185 %!error <Invalid call> sprandsym ()
14363
f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks
Rik <octave@nomad.inbox5.com>
parents: 14138
diff changeset
186 %!error sprandsym (ones (3), 0.5)
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
187 %!error sprandsym (3.5, 0.5)
28059
e9d57f6d6353 allow sprand and sprandsym to create empty sparse matrices
John W. Eaton <jwe@octave.org>
parents: 27923
diff changeset
188 %!error sprandsym (-1, 0.5)
13064
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
189 %!error sprandsym (3, -1)
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m
Rik <octave@nomad.inbox5.com>
parents: 11587
diff changeset
190 %!error sprandsym (3, 2)