octave-nkf: scripts/sparse/sprandsym.m annotate

annotate scripts/sparse/sprandsym.m @ 16036:dd583c258cf9

Fix trivial typo in run-octave.in comment

author	Jordi Gutiérrez Hermoso <jordigh@octave.org>
date	Mon, 11 Feb 2013 11:30:14 -0500
parents	86854d032a37
children	095d8d22c87d

rev	line source
14138 72c96de7a403 maint: update copyright notices for 2012 John W. Eaton <jwe@octave.org> parents: 13857 diff changeset	1 ## Copyright (C) 2004-2012 David Bateman and Andy Adler
72c96de7a403 maint: update copyright notices for 2012 John W. Eaton <jwe@octave.org> parents: 13857 diff changeset	2 ## Copyright (C) 2012 Jordi Gutiérrez Hermoso
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	3 ##
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6498 diff changeset	4 ## This file is part of Octave.
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	5 ##
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6498 diff changeset	6 ## Octave is free software; you can redistribute it and/or modify it
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6498 diff changeset	7 ## under the terms of the GNU General Public License as published by
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6498 diff changeset	8 ## the Free Software Foundation; either version 3 of the License, or (at
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6498 diff changeset	9 ## your option) any later version.
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6498 diff changeset	10 ##
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6498 diff changeset	11 ## 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	12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6498 diff changeset	13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6498 diff changeset	14 ## General Public License for more details.
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	15 ##
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	16 ## 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	17 ## along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6498 diff changeset	18 ## <http://www.gnu.org/licenses/>.
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 ## -- texinfo --
10793 be55736a0783 Grammarcheck the documentation from m-files. Rik <octave@nomad.inbox5.com> parents: 10549 diff changeset	21 ## @deftypefn {Function File} {} sprandsym (@var{n}, @var{d})
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	22 ## @deftypefnx {Function File} {} sprandsym (@var{s})
9051 1bf0ce0930be Grammar check TexInfo in all .m files Rik <rdrider0-list@yahoo.com> parents: 8920 diff changeset	23 ## Generate a symmetric random sparse matrix. The size of the matrix will be
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	24 ## @var{n} by @var{n}, with a density of values given by @var{d}.
14373 d00900b3dc4b doc: Use two spaces at start of sentence. Rik <octave@nomad.inbox5.com> parents: 14363 diff changeset	25 ## @var{d} should be between 0 and 1. Values will be normally
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	26 ## distributed with mean of zero and variance 1.
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	27 ##
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	28 ## If called with a single matrix argument, a random sparse matrix is
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	29 ## generated wherever the matrix @var{S} is non-zero in its lower
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	30 ## triangular part.
5642 2618a0750ae6 [project @ 2006-03-06 21:26:48 by jwe] jwe parents: 5610 diff changeset	31 ## @seealso{sprand, sprandn}
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	32 ## @end deftypefn
9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	33
6252 738c97e101eb [project @ 2007-01-24 19:58:46 by jwe] jwe parents: 6046 diff changeset	34 function S = sprandsym (n, d)
8506 bc982528de11 comment style fixes John W. Eaton <jwe@octave.org> parents: 7505 diff changeset	35
13064 bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	36 if (nargin != 1 && nargin != 2)
6046 34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe] jwe parents: 5642 diff changeset	37 print_usage ();
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	38 endif
13064 bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	39
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	40 if (nargin == 1)
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	41 [i, j] = find (tril (n));
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	42 [nr, nc] = size (n);
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	43 S = sparse (i, j, randn (size (i)), nr, nc);
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	44 S = S + tril (S, -1)';
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	45 return;
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	46 endif
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	47
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	48 if (!(isscalar (n) && n == fix (n) && n > 0))
13240 32980cbf2338 Use correct function name in error message in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13220 diff changeset	49 error ("sprandsym: N must be an integer greater than 0");
13064 bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	50 endif
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	51
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	52 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	53 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	54 endif
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	55
13205 abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	56 ## Actual number of nonzero entries
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	57 k = round (n^2*d);
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	58
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	59 ## Diagonal nonzero entries, same parity as k
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	60 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	61
13205 abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	62 ## Off diagonal nonzero entries
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	63 m = (k - r)/2;
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	64
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	65 ondiag = randperm (n, r);
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	66 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	67
13212 e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	68 ## Row index
13213 544304a09e42 Fix offbyones and typos in sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13212 diff changeset	69 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	70
13205 abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	71 ## Column index
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	72 j = offdiag - (i - 1).*(i - 2)/2;
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	73
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	74 diagvals = randn (1, r);
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	75 offdiagvals = randn (1, m);
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	76
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	77 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	78 [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	79
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	80 endfunction
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 function r = pick_rand_diag (n, k)
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	83 ## 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	84 ## 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	85 ## 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	86 ##
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	87 ## Let D be the number of diagonal entries and M the number of
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	88 ## off-diagonal entries. Then K = D + 2M. Let A = N(N-1)/2 be the
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	89 ## 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	90 ## 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	91 ##
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	92 ## T = nchoosek (N, D) * nchoosek (A, M)
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	93 ##
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	94 ## symmetric NxN matrices with a total of K nonzero entries and D on
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	95 ## the diagonal. Letting D range from mod (K,2) through min (N,K), and
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	96 ## 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	97 ## of those values.
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 ## However, we cannot use this form for computation, as the binomial
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	100 ## coefficients become unmanageably large. Instead, we use the
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	101 ## 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	102 ## be
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 ## (N - D)(N - D - 1)M
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	105 ## Q = -------------------------------
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	106 ## (D + 2)(D + 1)(A - M + 1)
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	107 ##
13212 e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	108 ## 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	109 ##
13212 e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	110 ## 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	111 ##
13212 e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	112 ## Their sum is thus S = sum (T)/T(1), and then C(i)/S is the desired
e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	113 ## probability P(i) for i=1:N. The cumsum will finally give the
e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	114 ## 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	115 ## the diagonal R.
13220 410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	116 ##
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	117 ## 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	118 ## 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	119 ## numerically to avoid overflow.
13212 e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	120
e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	121 ## Degenerate case
14552 86854d032a37 maint: miscellaneous style fixes for .m files John W. Eaton <jwe@octave.org> parents: 14373 diff changeset	122 if (k == 1)
13212 e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	123 r = 1;
e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	124 return
e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	125 endif
13064 bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	126
13205 abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	127 ## Compute the stuff described above
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	128 a = n*(n - 1)/2;
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	129 d = [mod(k,2):2:min(n,k)-2];
16036 dd583c258cf9 Fix trivial typo in run-octave.in comment Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14552 diff changeset	130
dd583c258cf9 Fix trivial typo in run-octave.in comment Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14552 diff changeset	131 ## Since we
dd583c258cf9 Fix trivial typo in run-octave.in comment Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14552 diff changeset	132 if (isempty (d))
dd583c258cf9 Fix trivial typo in run-octave.in comment Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14552 diff changeset	133 keyboard
dd583c258cf9 Fix trivial typo in run-octave.in comment Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14552 diff changeset	134 endif
dd583c258cf9 Fix trivial typo in run-octave.in comment Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 14552 diff changeset	135
13205 abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	136 m = (k - d)/2;
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	137 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	138
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	139 ## 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	140 ## 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	141 ## zero probabilities).
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	142 [~, midx] = max (cumsum (log (q))) ;
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	143 midx++;
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	144 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	145 rc = cumprod (q(midx:end));
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	146
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	147 ## 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	148 c = [lc, 1, rc];
13212 e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	149 s = sum (c);
13205 abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	150 p = c/s;
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	151
13212 e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	152 ## Add final d
13205 abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	153 d(end+1) = d(end) + 2;
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	154
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	155 ## Pick a random r using this distribution
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	156 r = d(sum (cumsum (p) < rand) + 1);
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	157
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	158 endfunction
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	159
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	160
13064 bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	161 %!test
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	162 %! s = sprandsym (10, 0.1);
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	163 %! assert (issparse (s));
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	164 %! assert (issymmetric (s));
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	165 %! assert (size (s), [10, 10]);
13205 abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	166 %! 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	167
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	168 %% Test 1-input calling form
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	169 %!test
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	170 %! 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	171 %! [i, j] = find (s);
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	172 %! assert (sort (i), [2 3]');
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	173 %! assert (sort (j), [2 3]');
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	174
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	175 %% Test input validation
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	176 %!error sprandsym ()
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	177 %!error sprandsym (1, 2, 3)
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	178 %!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	179 %!error sprandsym (3.5, 0.5)
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	180 %!error sprandsym (0, 0.5)
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	181 %!error sprandsym (3, -1)
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	182 %!error sprandsym (3, 2)
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	183

Mercurial > octave-nkf

annotate scripts/sparse/sprandsym.m @ 16036:dd583c258cf9