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

annotate scripts/sparse/sprandsym.m @ 20815:a260a6acb70f

fix test failures introduced by a22d8a2eb0e5 * scripts/ode/private/integrate_adaptive.m: fix stepping backwards, fix invocation of OutputFcn, fix text of some error messages * scripts/ode/private/integrate_const.m: remove use of option OutputSave * scripts/ode/private/integrate_n_steps.m: remove use of option OutputSave

author	Carlo de Falco <carlo.defalco@polimi.it>
date	Sun, 11 Oct 2015 23:09:01 +0200
parents	83792dd9bcc1
children

rev	line source
19898 4197fc428c7d maint: Update copyright notices for 2015. John W. Eaton <jwe@octave.org> parents: 19007 diff changeset	1 ## Copyright (C) 2004-2015 David Bateman and Andy Adler
14138 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})
19004 53af80da6781 doc: Update documentation of sparse functions including seealso links. Rik <rik@octave.org> parents: 17744 diff changeset	23 ## Generate a symmetric random sparse matrix.
5610 9761b7d24e9e [project @ 2006-02-09 09:12:02 by dbateman] dbateman parents: diff changeset	24 ##
19004 53af80da6781 doc: Update documentation of sparse functions including seealso links. Rik <rik@octave.org> parents: 17744 diff changeset	25 ## 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	26 ## 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	27 ## 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	28 ##
53af80da6781 doc: Update documentation of sparse functions including seealso links. Rik <rik@octave.org> parents: 17744 diff changeset	29 ## If called with a single matrix argument, a random sparse matrix is generated
20374 df437a52bcaf doc: Update more docstrings to have one sentence summary as first line. Rik <rik@octave.org> parents: 20038 diff changeset	30 ## wherever the matrix @var{s} is nonzero in its lower triangular part.
19004 53af80da6781 doc: Update documentation of sparse functions including seealso links. Rik <rik@octave.org> parents: 17744 diff changeset	31 ## @seealso{sprand, sprandn, spones, sparse}
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);
20441 83792dd9bcc1 Use in-place operators in m-files where possible. Rik <rik@octave.org> parents: 20374 diff changeset	44 S += tril (S, -1)';
13064 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;
17312 088d014a7fe2 Use semicolon after "return" statement in core m-files. Rik <rik@octave.org> parents: 16037 diff changeset	124 return;
13212 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];
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	130 m = (k - d)/2;
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	131 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	132
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	133 ## 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	134 ## 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	135 ## zero probabilities).
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	136 [~, midx] = max (cumsum (log (q))) ;
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	137 midx++;
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	138 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	139 rc = cumprod (q(midx:end));
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	140
410de573089a Avoid overflow in sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13213 diff changeset	141 ## 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	142 c = [lc, 1, rc];
13212 e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	143 s = sum (c);
13205 abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	144 p = c/s;
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	145
13212 e81b93284605 Various minor stylistic improvements to sprandsym.m Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13208 diff changeset	146 ## Add final d
13205 abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	147 d(end+1) = d(end) + 2;
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	148
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	149 ## Pick a random r using this distribution
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	150 r = d(sum (cumsum (p) < rand) + 1);
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	151
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	152 endfunction
abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	153
14363 f3d52523cde1 Use Octave coding conventions in all m-file %!test blocks Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	154
13064 bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	155 %!test
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	156 %! s = sprandsym (10, 0.1);
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	157 %! assert (issparse (s));
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	158 %! assert (issymmetric (s));
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	159 %! assert (size (s), [10, 10]);
13205 abf1e00111dd Completely new implementation of sprandsym Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 13064 diff changeset	160 %! 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	161
20038 9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions. Rik <rik@octave.org> parents: 19898 diff changeset	162 ## 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	163 %!test
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	164 %! 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	165 %! [i, j] = find (s);
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	166 %! assert (sort (i), [2 3]');
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	167 %! assert (sort (j), [2 3]');
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	168
20038 9fc020886ae9 maint: Clean up m-files to follow Octave coding conventions. Rik <rik@octave.org> parents: 19898 diff changeset	169 ## Test input validation
13064 bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	170 %!error sprandsym ()
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	171 %!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	172 %!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	173 %!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	174 %!error sprandsym (0, 0.5)
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	175 %!error sprandsym (3, -1)
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	176 %!error sprandsym (3, 2)
bae887ebea48 codesprint: Add input validation and tests for sprandsym.m Rik <octave@nomad.inbox5.com> parents: 11587 diff changeset	177

Mercurial > hg > octave-nkf

annotate scripts/sparse/sprandsym.m @ 20815:a260a6acb70f