[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Importing FIXES functions from octave-forge, and a couple of other thing
|
From: |
David Bateman |
|
Subject: |
Importing FIXES functions from octave-forge, and a couple of other things |
|
Date: |
Sun, 02 Apr 2006 23:37:33 +0200 |
|
User-agent: |
Mozilla Thunderbird 1.0.6-7.5.20060mdk (X11/20050322) |
I really don't like the presence of the FIXES directory in octave-forge
and would like to see it gone. In a previous mail on the octave-forge
list I discussed what I thought needed to be done to get the changed
bits in the octave-forge FIXES directory into octave, and attached here
is a first version of what I've come up with to do this.
The first patch is for scripts/plot/grid.m, where the octave-forge
version probes the state of the grid and toggles the grid in a similar
manner to matlab. This uses the octave-forge gget function to probe the
gnuplot grid state through a temporary file, which has potential race
issues. What I propose for grid.m is to use a persistent variable store
the state of the grid assuming it is initially off, and toggle in this
manner. The octave-forge version also allows minor grid lines, but with
a syntax that is different than matlab. I propose adding the "minor"
argument, but equally allow the following extensions to matlab
grid minor on
grid minor off
grid minor #
where # is the number of minor grid lines. Is it really necessary to
allow different number of minor grid lines in the different directions?
I didn't think so and so didn't allow it.
The second patch is to port the octave-forge random number generators to
octave, while keeping the old generators in place. The old generators
are selected by setting the "seed" of the generator with the "seed"
keyword, while the new generators are selected with the "state" keyword
and a vector with up to 625 elements. In addition the "twister" keyword
is an alais for the "state" keyword as R14sp3 allows this to select the
Mersenne Twister generator.
Additionally I included the rande, randg and randp generators and
include randlib versions of each of this as well. To test the code, the
test suite from the octave-forge code was also added. With this change
both the existing octave and octave-forge generators are available with
exactly the same sequences as previously. I good further suggestion is
that the statistics functions are updated to use the rande, randp and
randg functions to significantly accelerate their performance.
The major difference here with the old generators at the moment is that
the new generators all rely on a single mersenne twister, while teh old
generators each have a different generator. John what do you prefer.
The third patch adds an Fresize function. There were certain
difficulties with this. Firstly, the resize method of octave_value
didn't say whether the resize value was to be filled or not. The result
is that there are uninitialized values in the resize matrix. Therefore,
I modified the resize method of octave_value and all child sub-classes
to include a boolean argument "fill". Another problem was that the
resize function of Sparse<T> used ridx et al, rather than xridx and so
make_unique was called, when the old rep save already saved. The result
was that the rep->count on the original matrix was decremented twice and
the original matrix was deleted!!!
Now the reason for the Fresize function was to allow
retval = zeros (nr, nc, class (x));
in triu, tril, hankel.m and toeplitz functions in FIXES to be replaced with
retval = resize (resize (x, 0), nr, nc);
so that these functions would work correctly with user types in a
similar manner to the octave-forge functions but without the speed
memory trade-off made by the octave-forge functions. That is tril, etc
can't work with an arbitrary type as written in octave.
The fourth patch was as recently discussed I added function handles and
cell arrays of function handles to the fsolve, dasrt, etc functions.
This allows fixed arguments to be handled like "fsolve (@(x) myfun
(x,5), x_init)". Another small change I made to these functions is that
is the function was previously given a string it created a function in
the symbol table called fsolve_fcn and perhaps also fsolve_jac. However,
these were never deleted. This has two problems. Firstly the population
of the symbol table, and secondly if we ever compile the underlying
fortran code to allow recursion, then this was non reentrant. So I
created a unique function name for each function and deleted them when
leaving fsolve.
A fifth minor patch was to get rid of the compile warning in ls-hdf5.h,
and reverse the sense of the isscalar and isvector tests in test_number.m.
The final patch included here is the fftw-wisdom patch I sent last week
as it doesn't appear to be in the CVS (updated with John's comment).
All of these patches should allow all but tf2zp.m and zp2tf.m to be
deleted from the octave-forge FIXES directory if they are accepted. As
for zp2tf.m and tf2zp.m I still don't know what bugs these versions are
supposedly fixing.. The only thing I can see related to this is
http://velveeta.che.wisc.edu/octave/lists/archive//bug-octave.2001/msg00171.html
and the error in this message now doesn't appear to happen with the
octave versions. So perhaps Paul you can confirm, but I think zp2tf.m
and tf2zp.m can also be deleted. John, comments please and which patches
need more work before you'll accept them and what is that work?
Cheers, enjoy.
David
--
David Bateman address@hidden
Motorola Labs - Paris +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin +33 6 72 01 06 33 (Mob)
91193 Gif-Sur-Yvette FRANCE +33 1 69 35 77 01 (Fax)
The information contained in this communication has been classified as:
[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary
2006-04-02 David Bateman <address@hidden>
* plot/grid.m: cache the state of the grid to allow
toggling. Add keyword "minor" for minor grid.
2006-04-02 David Bateman <address@hidden>
* DLD-FUNCTIONS/daspk.cc (Fdaspk), DLD-FUNCTIONS/daspk.cc
(Fdasrt), DLD-FUNCTIONS/daspk.cc (Fdassl),
DLD-FUNCTIONS/daspk.cc (Ffsolve), DLD-FUNCTIONS/daspk.cc
(Flsode): Add the ability to read function handles, inline
functions, and cell arrays of strings, inline and function
handles.
2006-04-02 David Bateman <address@hidden>
* ls-hdf5.h (open): remove unused arg prot.
2006-04-02 David Bateman <address@hidden>
* test_number.m: Reverse sense of isscalar and isvector tests
for recent changes.
2006-04-02 David Bateman <address@hidden>
* interpreter/matrix.txi: Add rande, randp, randg and update
for different random generator behavior.
2006-04-02 David Bateman <address@hidden>
* randlib/wrap.f: Add genexp, gengam and ignpoi functions to
the wrapped functions.
2006-04-02 David Bateman <address@hidden>
* Makefile.in (INCLUDES): Add randgamma.h, randpoisson.h and
randmtzig.h.
(LIBOCTAVE_C_SOURCES): Add randgamma.c, randpoisson.c and
randmtzig.c.
* oct-rand.cc (do_new_initialization): Rename from
do_initialization.
(maybe_initialize): Initialize new generators as well.
(ColumnVectoroctave_rand::state (void), void
octave_rand::state (const ColumnVector &)): New functions.
(std::string octave_rand::distribution (void), void
octave_rand::distribution (const std::string &)): Add gamma,
exponential and poisson distributions.
(void octave_rand::exponential_distribution,
void octave_rand::poisson_distribution,
void octave_rand::gamma_distribution): New methods to select
exponential, poisson or gamma distribution.
(double octave_rand::scalar (void)), Matrix
octave_rand::matrix (octave_idx_type, octave_idx_type),
NDArray octave_rand::nd_array (const dim_vector &): Add new
distributions.
* oct-rand.h (static ColumnVector state (void), static void
state (const ColumnVector &), static void
exponential_distribution (void), poisson_distribution (void),
static void gamma_distribution (void)): New functions.
(static Matrix matrix (..., double), static NDArray nd_array
(..., double), static double scalar (double), static
Array<double> vector (..., double)): Additional argument for
gamma and poisson distributions.
* randpoisson.c, randpoisson.h, randgamma.c, randmtzig.c,
randmtzig.h: New files.
2006-04-02 David Bateman <address@hidden>
* DLD-FUNCTIONS/rand.cc (do_rand): Additional argument for
gamma and poisson distributions. Change "state" and "seed"
arguments so that they choose between generators.
Add, poisson, gamma and exponential generators.
(Frand Frandn): Update docs for new generators, add tests.
(Frande, Frandp, Frandg): New generators, with test code.
2006-04-02 David Bateman <address@hidden>
* Sparse.cc (resize): Use xcidx rather than cdix, etc to avoid
copy of original matrix.
2006-04-02 David Bateman <address@hidden>
* general/tril.m, general.triu.m: Use resize(resize(x, 0),
nr, nc) rather than zeros(nr, nc) to allow user types to
work correctly.
* special-matrix/hankel.m, special-matrix/toeplitz.m: ditto.
2006-04-02 David Bateman <address@hidden>
* data.cc (Fresize): New function.
* oct-map.cc, ov-base-mat.cc, ov-base-sparse.cc, ov-base.cc,
ov-bool.cc, ov-complex.cc, ov-range.cc, ov-scalar.cc,
ov-str-mat.cc (resize): Add boolean fill argument.
* oct-map.h, ov-base-mat.h, ov-base-sparse.h, ov-base.h,
ov-bool.h, ov-complex.h, ov-intx.h, ov-range.h, ov-scalar.h,
ov-str-mat.h, ov.h (resize): ditto.
2006-03-30 David Bateman <address@hidden>
* DLD-FUNCTIONS/fftw_wisdom.cc: Don't attempt to save wisdom to
an empty filename or invalid filename.
Index: scripts/plot/grid.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/plot/grid.m,v
retrieving revision 1.29
diff -c -r1.29 grid.m
*** scripts/plot/grid.m 6 Mar 2006 21:26:50 -0000 1.29
--- scripts/plot/grid.m 2 Apr 2006 20:23:27 -0000
***************
*** 19,27 ****
## -*- texinfo -*-
## @deftypefn {Function File} {} grid (@var{arg})
## For two-dimensional plotting, force the display of a grid on the plot.
## The argument may be either @code{"on"} or @code{"off"}. If it is
! ## omitted, @code{"on"} is assumed.
## @seealso{plot, semilogx, semilogy, loglog, polar, mesh, contour,
## bar, stairs, replot, xlabel, ylabel, title}
## @end deftypefn
--- 19,34 ----
## -*- texinfo -*-
## @deftypefn {Function File} {} grid (@var{arg})
+ ## @deftypefnx {Function File} {} grid ("minor", @var{arg2})
## For two-dimensional plotting, force the display of a grid on the plot.
## The argument may be either @code{"on"} or @code{"off"}. If it is
! ## omitted, the the current grid state is toggled.
! ##
! ## If @var{arg} is @code{"minor"} then the minor grid is toggled. When
! ## using a minor grid a second argument @var{arg2} is allowed, which can
! ## be either @code{"on"} or @code{"off"} to explicitly set the state of
! ## the minor grid, or alternatively a positive integer specifying the
! ## number of minor grid lines.
## @seealso{plot, semilogx, semilogy, loglog, polar, mesh, contour,
## bar, stairs, replot, xlabel, ylabel, title}
## @end deftypefn
***************
*** 30,51 ****
## PKG_ADD: mark_as_command grid
! function grid (x)
!
usage_msg = "grid (\"on\" | \"off\")";
do_replot = false;
if (nargin == 0)
! __gnuplot_raw__ ("set grid;\n");
do_replot = true;
elseif (nargin == 1)
if (ischar (x))
if (strcmp ("off", x))
__gnuplot_raw__ ("set nogrid;\n");
do_replot = true;
elseif (strcmp ("on", x))
__gnuplot_raw__ ("set grid;\n");
do_replot = true;
else
usage (usage_msg);
--- 37,83 ----
## PKG_ADD: mark_as_command grid
! function grid (x,y)
! persistent grid_on = false;
! persistent minor_on = false;
! persistent minor_tics = 5;
usage_msg = "grid (\"on\" | \"off\")";
do_replot = false;
if (nargin == 0)
! if ((grid_on = !grid_on))
! __gnuplot_raw__ ("set grid;\n");
! else
! __gnuplot_raw__ ("set nogrid;\n");
! endif
do_replot = true;
elseif (nargin == 1)
if (ischar (x))
if (strcmp ("off", x))
__gnuplot_raw__ ("set nogrid;\n");
+ grid_on = false;
do_replot = true;
elseif (strcmp ("on", x))
__gnuplot_raw__ ("set grid;\n");
+ grid_on = true;
+ do_replot = true;
+ elseif (strcmp ("minor", x))
+ if ((minor_on = !minor_on))
+ cmd = sprintf("set mxtics %d;", minor_tics);
+ __gnuplot_raw__ (cmd);
+ cmd = sprintf("set mytics %d;", minor_tics);
+ __gnuplot_raw__ (cmd);
+ __gnuplot_raw__ ("set grid xtics mxtics ytics mxtics;\n");
+ minor_on = true;
+ else
+ if (grid_on)
+ __gnuplot_raw__ ("set grid xtics nomxtics ytics nomxtics;\n");
+ else
+ __gnuplot_raw__ ("set grid noxtics nomxtics noytics nomxtics;\n");
+ endif
+ minor_on = false;
+ endif
do_replot = true;
else
usage (usage_msg);
***************
*** 53,61 ****
else
error ("grid: argument must be a string");
endif
else
usage (usage_msg);
! endif
if (do_replot && automatic_replot)
replot ();
--- 85,136 ----
else
error ("grid: argument must be a string");
endif
+ elseif (nargin == 2)
+ if (ischar (x))
+ if (strcmp ("minor", x))
+ d = str2num (y);
+ if (isempty(d))
+ if (strcmp ("off", y))
+ if (grid_on)
+ __gnuplot_raw__ ("set grid xtics nomxtics ytics nomxtics;\n");
+ else
+ __gnuplot_raw__ ("set grid noxtics nomxtics noytics nomxtics;\n");
+ endif
+ minor_on = false;
+ elseif (strcmp ("on", y))
+ cmd = sprintf("set mxtics %d;", minor_tics);
+ __gnuplot_raw__ (cmd);
+ cmd = sprintf("set mytics %d;", minor_tics);
+ __gnuplot_raw__ (cmd);
+ __gnuplot_raw__ ("set grid xtics mxtics ytics mxtics;\n");
+ minor_on = true;
+ else
+ usage (usage_msg);
+ endif
+ do_replot = true;
+ else
+ if (isscalar(d) && !isnan(d) && !isinf(d))
+ minor_tics = max(floor(d), 0);
+ cmd = sprintf("set mxtics %d;", minor_tics);
+ __gnuplot_raw__ (cmd);
+ cmd = sprintf("set mytics %d;", minor_tics);
+ __gnuplot_raw__ (cmd);
+ __gnuplot_raw__ ("set grid xtics mxtics ytics mxtics;\n");
+ minor_on = true;
+ do_replot = true;
+ else
+ usage (usage_msg);
+ endif
+ endif
+ else
+ usage (usage_msg);
+ endif
+ else
+ usage (usage_msg);
+ endif
else
usage (usage_msg);
! endif
if (do_replot && automatic_replot)
replot ();
Index: src/DLD-FUNCTIONS/daspk.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/DLD-FUNCTIONS/daspk.cc,v
retrieving revision 1.15
diff -c -r1.15 daspk.cc
*** src/DLD-FUNCTIONS/daspk.cc 6 Mar 2006 21:26:54 -0000 1.15
--- src/DLD-FUNCTIONS/daspk.cc 2 Apr 2006 20:23:31 -0000
***************
*** 37,42 ****
--- 37,43 ----
#include "gripes.h"
#include "oct-obj.h"
#include "ov-fcn.h"
+ #include "ov-cell.h"
#include "pager.h"
#include "unwind-prot.h"
#include "utils.h"
***************
*** 205,213 ****
row of the output @var{x} is @var{x_0} and the first row\n\
of the output @var{xdot} is @var{xdot_0}.\n\
\n\
! The first argument, @var{fcn}, is a string that names the function to\n\
! call to compute the vector of residuals for the set of equations.\n\
! It must have the form\n\
\n\
@example\n\
@var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
--- 206,215 ----
row of the output @var{x} is @var{x_0} and the first row\n\
of the output @var{xdot} is @var{xdot_0}.\n\
\n\
! The first argument, @var{fcn}, is a string or a two element cell array\n\
! of strings, inline or function handle, that names the function, to call\n\
! to compute the vector of residuals for the set of equations. It must\n\
! have the form\n\
\n\
@example\n\
@var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
***************
*** 286,328 ****
if (nargin > 3 && nargin < 6)
{
daspk_fcn = 0;
daspk_jac = 0;
octave_value f_arg = args(0);
! switch (f_arg.rows ())
{
! case 1:
! daspk_fcn = extract_function
! (args(0), "daspk", "__daspk_fcn__",
! "function res = __daspk_fcn__ (x, xdot, t) res = ",
! "; endfunction");
! break;
!
! case 2:
! {
! string_vector tmp = f_arg.all_strings ();
!
! if (! error_state)
! {
! daspk_fcn = extract_function
! (tmp(0), "daspk", "__daspk_fcn__",
! "function res = __daspk_fcn__ (x, xdot, t) res = ",
! "; endfunction");
! if (daspk_fcn)
{
! daspk_jac = extract_function
! (tmp(1), "daspk", "__daspk_jac__",
! "function jac = __daspk_jac__ (x, xdot, t, cj) jac = ",
! "; endfunction");
! if (! daspk_jac)
! daspk_fcn = 0;
}
! }
! }
}
if (error_state || ! daspk_fcn)
--- 288,399 ----
if (nargin > 3 && nargin < 6)
{
+ std::string fcn_name, fname, jac_name, jname;
daspk_fcn = 0;
daspk_jac = 0;
octave_value f_arg = args(0);
! if (f_arg.is_cell ())
! {
! Cell c = f_arg.cell_value ();
! if (c.length() == 1)
! f_arg = c(0);
! else if (c.length() == 2)
! {
! if (c(0).is_function_handle () || c(0).is_inline_function ())
! daspk_fcn = c(0).function_value ();
! else
! {
! fcn_name = unique_symbol_name ("__daspk_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, xdot, t) y = ");
! daspk_fcn = extract_function
! (c(0), "daspk", fcn_name, fname, "; endfunction");
! }
!
! if (daspk_fcn)
! {
! if (c(1).is_function_handle () || c(1).is_inline_function ())
! daspk_jac = c(1).function_value ();
! else
! {
! jac_name = unique_symbol_name ("__daspk_jac__");
! jname = "function jac = ";
! jname.append(jac_name);
! jname.append (" (x, xdot, t, cj) jac = ");
! daspk_jac = extract_function
! (c(1), "daspk", jac_name, jname, "; endfunction");
!
! if (!daspk_jac)
! {
! if (fcn_name.length())
! clear_function (fcn_name);
! daspk_fcn = 0;
! }
! }
! }
! }
! else
! DASPK_ABORT1 ("incorrect number of elements in cell array");
! }
!
! if (!daspk_fcn && ! f_arg.is_cell())
{
! if (f_arg.is_function_handle () || f_arg.is_inline_function ())
! daspk_fcn = f_arg.function_value ();
! else
! {
! switch (f_arg.rows ())
! {
! case 1:
! do
! {
! fcn_name = unique_symbol_name ("__daspk_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, xdot, t) y = ");
! daspk_fcn = extract_function
! (f_arg, "daspk", fcn_name, fname, "; endfunction");
! }
! while (0);
! break;
! case 2:
{
! string_vector tmp = f_arg.all_strings ();
! if (! error_state)
! {
! fcn_name = unique_symbol_name ("__daspk_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, xdot, t) y = ");
! daspk_fcn = extract_function
! (tmp(0), "daspk", fcn_name, fname, "; endfunction");
!
! if (daspk_fcn)
! {
! jac_name = unique_symbol_name ("__daspk_jac__");
! jname = "function jac = ";
! jname.append(jac_name);
! jname.append (" (x, xdot, t, cj) jac = ");
! daspk_jac = extract_function
! (tmp(1), "daspk", jac_name, jname,
! "; endfunction");
!
! if (!daspk_jac)
! {
! if (fcn_name.length())
! clear_function (fcn_name);
! daspk_fcn = 0;
! }
! }
! }
}
! }
! }
}
if (error_state || ! daspk_fcn)
***************
*** 375,380 ****
--- 446,456 ----
else
output = dae.integrate (out_times, deriv_output);
+ if (fcn_name.length())
+ clear_function (fcn_name);
+ if (jac_name.length())
+ clear_function (jac_name);
+
if (! error_state)
{
std::string msg = dae.error_message ();
Index: src/DLD-FUNCTIONS/dasrt.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/DLD-FUNCTIONS/dasrt.cc,v
retrieving revision 1.19
diff -c -r1.19 dasrt.cc
*** src/DLD-FUNCTIONS/dasrt.cc 20 Mar 2006 18:29:20 -0000 1.19
--- src/DLD-FUNCTIONS/dasrt.cc 2 Apr 2006 20:23:31 -0000
***************
*** 36,41 ****
--- 36,42 ----
#include "gripes.h"
#include "oct-obj.h"
#include "ov-fcn.h"
+ #include "ov-cell.h"
#include "pager.h"
#include "parse.h"
#include "unwind-prot.h"
***************
*** 249,257 ****
@var{t_out} will be the point at which the stopping condition was met,\n\
and may not correspond to any element of the vector @var{t}.\n\
\n\
! The first argument, @var{fcn}, is a string that names the function to\n\
! call to compute the vector of residuals for the set of equations.\n\
! It must have the form\n\
\n\
@example\n\
@var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
--- 250,258 ----
@var{t_out} will be the point at which the stopping condition was met,\n\
and may not correspond to any element of the vector @var{t}.\n\
\n\
! The first argument, @var{fcn}, is a string, or cell array of strings or\n\
! inline or function handles, that names the function to call to compute\n\
! the vector of residuals for the set of equations. It must have the form\n\
\n\
@example\n\
@var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
***************
*** 261,269 ****
in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\
scalar.\n\
\n\
! If @var{fcn} is a two-element string array, the first element names\n\
! the function @math{f} described above, and the second element names\n\
! a function to compute the modified Jacobian\n\
\n\
@tex\n\
$$\n\
--- 262,270 ----
in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\
scalar.\n\
\n\
! If @var{fcn} is a two-element string array, or two element cell array,\n\
! the first element names the function @math{f} described above, and the\n\
! second element names a function to compute the modified Jacobian\n\
\n\
@tex\n\
$$\n\
***************
*** 367,372 ****
--- 368,374 ----
return retval;
}
+ std::string fcn_name, fname, jac_name, jname;
dasrt_f = 0;
dasrt_j = 0;
dasrt_cf = 0;
***************
*** 377,419 ****
octave_value f_arg = args(0);
! switch (f_arg.rows ())
{
! case 1:
! dasrt_f = extract_function
! (args(0), "dasrt", "__dasrt_fcn__",
! "function res = __dasrt_fcn__ (x, xdot, t) res = ",
! "; endfunction");
! break;
! case 2:
! {
! string_vector tmp = args(0).all_strings ();
! if (! error_state)
! {
! dasrt_f = extract_function
! (tmp(0), "dasrt", "__dasrt_fcn__",
! "function res = __dasrt_fcn__ (x, xdot, t) res = ",
! "; endfunction");
! if (dasrt_f)
! {
! dasrt_j = extract_function
! (tmp(1), "dasrt", "__dasrt_jac__",
! "function jac = __dasrt_jac__ (x, xdot, t, cj) jac = ",
! "; endfunction");
!
! if (! dasrt_j)
! dasrt_f = 0;
}
! }
! }
! break;
! default:
! DASRT_ABORT1
! ("first arg should be a string or 2-element string array");
}
if (error_state || (! dasrt_f))
--- 379,480 ----
octave_value f_arg = args(0);
! if (f_arg.is_cell ())
{
! Cell c = f_arg.cell_value ();
! if (c.length() == 1)
! f_arg = c(0);
! else if (c.length() == 2)
! {
! if (c(0).is_function_handle () || c(0).is_inline_function ())
! dasrt_f = c(0).function_value ();
! else
! {
! fcn_name = unique_symbol_name ("__dasrt_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, xdot, t) y = ");
! dasrt_f = extract_function
! (c(0), "dasrt", fcn_name, fname, "; endfunction");
! }
!
! if (dasrt_f)
! {
! if (c(1).is_function_handle () || c(1).is_inline_function ())
! dasrt_j = c(1).function_value ();
! else
! {
! jac_name = unique_symbol_name ("__dasrt_jac__");
! jname = "function jac = ";
! jname.append(jac_name);
! jname.append (" (x, xdot, t, cj) jac = ");
! dasrt_j = extract_function
! (c(1), "dasrt", jac_name, jname, "; endfunction");
!
! if (!dasrt_j)
! {
! if (fcn_name.length())
! clear_function (fcn_name);
! dasrt_f = 0;
! }
! }
! }
! }
! else
! DASRT_ABORT1 ("incorrect number of elements in cell array");
! }
!
! if (!dasrt_f && ! f_arg.is_cell())
! {
! if (f_arg.is_function_handle () || f_arg.is_inline_function ())
! dasrt_f = f_arg.function_value ();
! else
! {
! switch (f_arg.rows ())
! {
! case 1:
! fcn_name = unique_symbol_name ("__dasrt_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, xdot, t) y = ");
! dasrt_f = extract_function
! (f_arg, "dasrt", fcn_name, fname, "; endfunction");
! break;
! case 2:
! {
! string_vector tmp = args(0).all_strings ();
! if (! error_state)
! {
! fcn_name = unique_symbol_name ("__dasrt_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, xdot, t) y = ");
! dasrt_f = extract_function
! (tmp(0), "dasrt", fcn_name, fname, "; endfunction");
! if (dasrt_f)
! {
! jac_name = unique_symbol_name ("__dasrt_jac__");
! jname = "function jac = ";
! jname.append(jac_name);
! jname.append (" (x, xdot, t, cj) jac = ");
! dasrt_j = extract_function
! (tmp(1), "dasrt", jac_name, jname, "; endfunction");
!
! if (! dasrt_j)
! dasrt_f = 0;
! }
! }
}
! break;
! default:
! DASRT_ABORT1
! ("first arg should be a string or 2-element string array");
! }
! }
}
if (error_state || (! dasrt_f))
***************
*** 423,432 ****
argp++;
! if (args(1).is_string ())
{
! dasrt_cf = is_valid_function (args(1), "dasrt", true);
if (! dasrt_cf)
DASRT_ABORT1 ("expecting function name as argument 2");
--- 484,503 ----
argp++;
! if (args(1).is_function_handle() || args(1).is_inline_function())
{
! dasrt_cf = args(1).function_value();
!
! if (! dasrt_cf)
! DASRT_ABORT1 ("expecting function name as argument 2");
+ argp++;
+
+ func.set_constraint_function (dasrt_user_cf);
+ }
+ else if (args(1).is_string ())
+ {
+ dasrt_cf = is_valid_function (args(1), "dasrt", true);
if (! dasrt_cf)
DASRT_ABORT1 ("expecting function name as argument 2");
***************
*** 483,488 ****
--- 554,564 ----
else
output = dae.integrate (out_times);
+ if (fcn_name.length())
+ clear_function (fcn_name);
+ if (jac_name.length())
+ clear_function (jac_name);
+
if (! error_state)
{
std::string msg = dae.error_message ();
Index: src/DLD-FUNCTIONS/dassl.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/DLD-FUNCTIONS/dassl.cc,v
retrieving revision 1.33
diff -c -r1.33 dassl.cc
*** src/DLD-FUNCTIONS/dassl.cc 20 Mar 2006 18:30:21 -0000 1.33
--- src/DLD-FUNCTIONS/dassl.cc 2 Apr 2006 20:23:31 -0000
***************
*** 37,42 ****
--- 37,43 ----
#include "gripes.h"
#include "oct-obj.h"
#include "ov-fcn.h"
+ #include "ov-cell.h"
#include "pager.h"
#include "unwind-prot.h"
#include "utils.h"
***************
*** 208,216 ****
row of the output @var{x} is @var{x_0} and the first row\n\
of the output @var{xdot} is @var{xdot_0}.\n\
\n\
! The first argument, @var{fcn}, is a string that names the function to\n\
! call to compute the vector of residuals for the set of equations.\n\
! It must have the form\n\
\n\
@example\n\
@var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
--- 209,218 ----
row of the output @var{x} is @var{x_0} and the first row\n\
of the output @var{xdot} is @var{xdot_0}.\n\
\n\
! The first argument, @var{fcn}, is a string or a two element cell array\n\
! of strings, inline or function handle, that names the function, to call\n\
! to compute the vector of residuals for the set of equations. It must\n\
! have the form\n\
\n\
@example\n\
@var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
***************
*** 291,333 ****
if (nargin > 3 && nargin < 6 && nargout < 5)
{
dassl_fcn = 0;
dassl_jac = 0;
octave_value f_arg = args(0);
! switch (f_arg.rows ())
{
! case 1:
! dassl_fcn = extract_function
! (f_arg, "dassl", "__dassl_fcn__",
! "function res = __dassl_fcn__ (x, xdot, t) res = ",
! "; endfunction");
! break;
!
! case 2:
! {
! string_vector tmp = f_arg.all_strings ();
!
! if (! error_state)
! {
! dassl_fcn = extract_function
! (tmp(0), "dassl", "__dassl_fcn__",
! "function res = __dassl_fcn__ (x, xdot, t) res = ",
! "; endfunction");
! if (dassl_fcn)
{
! dassl_jac = extract_function
! (tmp(1), "dassl", "__dassl_jac__",
! "function jac = __dassl_jac__ (x, xdot, t, cj) jac = ",
! "; endfunction");
! if (! dassl_jac)
! dassl_fcn = 0;
}
! }
! }
}
if (error_state || ! dassl_fcn)
--- 293,404 ----
if (nargin > 3 && nargin < 6 && nargout < 5)
{
+ std::string fcn_name, fname, jac_name, jname;
dassl_fcn = 0;
dassl_jac = 0;
octave_value f_arg = args(0);
! if (f_arg.is_cell ())
! {
! Cell c = f_arg.cell_value ();
! if (c.length() == 1)
! f_arg = c(0);
! else if (c.length() == 2)
! {
! if (c(0).is_function_handle () || c(0).is_inline_function ())
! dassl_fcn = c(0).function_value ();
! else
! {
! fcn_name = unique_symbol_name ("__dassl_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, xdot, t) y = ");
! dassl_fcn = extract_function
! (c(0), "dassl", fcn_name, fname, "; endfunction");
! }
!
! if (dassl_fcn)
! {
! if (c(1).is_function_handle () || c(1).is_inline_function ())
! dassl_jac = c(1).function_value ();
! else
! {
! jac_name = unique_symbol_name ("__dassl_jac__");
! jname = "function jac = ";
! jname.append(jac_name);
! jname.append (" (x, xdot, t, cj) jac = ");
! dassl_jac = extract_function
! (c(1), "dassl", jac_name, jname, "; endfunction");
!
! if (!dassl_jac)
! {
! if (fcn_name.length())
! clear_function (fcn_name);
! dassl_fcn = 0;
! }
! }
! }
! }
! else
! DASSL_ABORT1 ("incorrect number of elements in cell array");
! }
!
! if (!dassl_fcn && ! f_arg.is_cell())
{
! if (f_arg.is_function_handle () || f_arg.is_inline_function ())
! dassl_fcn = f_arg.function_value ();
! else
! {
! switch (f_arg.rows ())
! {
! case 1:
! do
! {
! fcn_name = unique_symbol_name ("__dassl_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, xdot, t) y = ");
! dassl_fcn = extract_function
! (f_arg, "dassl", fcn_name, fname, "; endfunction");
! }
! while (0);
! break;
! case 2:
{
! string_vector tmp = f_arg.all_strings ();
! if (! error_state)
! {
! fcn_name = unique_symbol_name ("__dassl_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, xdot, t) y = ");
! dassl_fcn = extract_function
! (tmp(0), "dassl", fcn_name, fname, "; endfunction");
!
! if (dassl_fcn)
! {
! jac_name = unique_symbol_name ("__dassl_jac__");
! jname = "function jac = ";
! jname.append(jac_name);
! jname.append (" (x, xdot, t, cj) jac = ");
! dassl_jac = extract_function
! (tmp(1), "dassl", jac_name, jname,
! "; endfunction");
!
! if (!dassl_jac)
! {
! if (fcn_name.length())
! clear_function (fcn_name);
! dassl_fcn = 0;
! }
! }
! }
}
! }
! }
}
if (error_state || ! dassl_fcn)
***************
*** 381,386 ****
--- 452,462 ----
else
output = dae.integrate (out_times, deriv_output);
+ if (fcn_name.length())
+ clear_function (fcn_name);
+ if (jac_name.length())
+ clear_function (jac_name);
+
if (! error_state)
{
std::string msg = dae.error_message ();
Index: src/DLD-FUNCTIONS/fsolve.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/DLD-FUNCTIONS/fsolve.cc,v
retrieving revision 1.27
diff -c -r1.27 fsolve.cc
*** src/DLD-FUNCTIONS/fsolve.cc 26 Apr 2005 19:24:35 -0000 1.27
--- src/DLD-FUNCTIONS/fsolve.cc 2 Apr 2006 20:23:31 -0000
***************
*** 37,42 ****
--- 37,43 ----
#include "gripes.h"
#include "oct-obj.h"
#include "ov-fcn.h"
+ #include "ov-cell.h"
#include "pager.h"
#include "unwind-prot.h"
#include "utils.h"
***************
*** 220,229 ****
and an initial starting point @var{x0}, @code{fsolve} solves the set of\n\
equations such that @code{f(@var{x}) == 0}.\n\
\n\
! If @var{fcn} is a two-element string array, the first element names\n\
! the function @math{f} described above, and the second element names\n\
! a function of the form @code{j (@var{x})} to compute the Jacobian\n\
! matrix with elements\n\
@tex\n\
$$ J = {\\partial f_i \\over \\partial x_j} $$\n\
@end tex\n\
--- 221,231 ----
and an initial starting point @var{x0}, @code{fsolve} solves the set of\n\
equations such that @code{f(@var{x}) == 0}.\n\
\n\
! If @var{fcn} is a two-element string array, or a two element cell array\n\
! containing either the function name or inline or function handle. The\n\
! first element names the function @math{f} described above, and the second\n\
! element names a function of the form @code{j (@var{x})} to compute the\n\
! Jacobian matrix with elements\n\
@tex\n\
$$ J = {\\partial f_i \\over \\partial x_j} $$\n\
@end tex\n\
***************
*** 257,299 ****
if (nargin == 2 && nargout < 4)
{
fsolve_fcn = 0;
fsolve_jac = 0;
octave_value f_arg = args(0);
! switch (f_arg.rows ())
{
! case 1:
! fsolve_fcn = extract_function
! (f_arg, "fsolve", "__fsolve_fcn__",
! "function y = __fsolve_fcn__ (x) y = ",
! "; endfunction");
! break;
!
! case 2:
! {
! string_vector tmp = f_arg.all_strings ();
!
! if (! error_state)
! {
! fsolve_fcn = extract_function
! (tmp(0), "fsolve", "__fsolve_fcn__",
! "function y = __fsolve_fcn__ (x) y = ",
! "; endfunction");
! if (fsolve_fcn)
{
! fsolve_jac = extract_function
! (tmp(1), "fsolve", "__fsolve_jac__",
! "function jac = __fsolve_jac__ (x) jac = ",
! "; endfunction");
! if (! fsolve_jac)
! fsolve_fcn = 0;
}
! }
! }
}
if (error_state || ! fsolve_fcn)
--- 259,370 ----
if (nargin == 2 && nargout < 4)
{
+ std::string fcn_name, fname, jac_name, jname;
fsolve_fcn = 0;
fsolve_jac = 0;
octave_value f_arg = args(0);
! if (f_arg.is_cell ())
! {
! Cell c = f_arg.cell_value ();
! if (c.length() == 1)
! f_arg = c(0);
! else if (c.length() == 2)
! {
! if (c(0).is_function_handle () || c(0).is_inline_function ())
! fsolve_fcn = c(0).function_value ();
! else
! {
! fcn_name = unique_symbol_name ("__fsolve_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x) y = ");
! fsolve_fcn = extract_function
! (c(0), "fsolve", fcn_name, fname, "; endfunction");
! }
!
! if (fsolve_fcn)
! {
! if (c(1).is_function_handle () || c(1).is_inline_function ())
! fsolve_jac = c(1).function_value ();
! else
! {
! jac_name = unique_symbol_name ("__fsolve_jac__");
! jname = "function y = ";
! jname.append (jac_name);
! jname.append (" (x) jac = ");
! fsolve_jac = extract_function
! (c(1), "fsolve", jac_name, jname, "; endfunction");
!
! if (!fsolve_jac)
! {
! if (fcn_name.length())
! clear_function (fcn_name);
! fsolve_fcn = 0;
! }
! }
! }
! }
! else
! FSOLVE_ABORT1 ("incorrect number of elements in cell array");
! }
!
! if (!fsolve_fcn && ! f_arg.is_cell())
{
! if (f_arg.is_function_handle () || f_arg.is_inline_function ())
! fsolve_fcn = f_arg.function_value ();
! else
! {
! switch (f_arg.rows ())
! {
! case 1:
! do
! {
! fcn_name = unique_symbol_name ("__fsolve_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x) y = ");
! fsolve_fcn = extract_function
! (f_arg, "fsolve", fcn_name, fname, "; endfunction");
! }
! while (0);
! break;
! case 2:
{
! string_vector tmp = f_arg.all_strings ();
! if (! error_state)
! {
! fcn_name = unique_symbol_name ("__fsolve_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x) y = ");
! fsolve_fcn = extract_function
! (tmp(0), "fsolve", fcn_name, fname, "; endfunction");
!
! if (fsolve_fcn)
! {
! jac_name = unique_symbol_name ("__fsolve_jac__");
! jname = "function y = ";
! jname.append (jac_name);
! jname.append (" (x) jac = ");
! fsolve_jac = extract_function
! (tmp(1), "fsolve", jac_name, jname,
! "; endfunction");
!
! if (!fsolve_jac)
! {
! if (fcn_name.length())
! clear_function (fcn_name);
! fsolve_fcn = 0;
! }
! }
! }
}
! }
! }
}
if (error_state || ! fsolve_fcn)
***************
*** 320,325 ****
--- 391,401 ----
octave_idx_type info;
ColumnVector soln = nleqn.solve (info);
+ if (fcn_name.length())
+ clear_function (fcn_name);
+ if (jac_name.length())
+ clear_function (jac_name);
+
if (! error_state)
{
std::string msg = nleqn.error_message ();
Index: src/DLD-FUNCTIONS/lsode.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/DLD-FUNCTIONS/lsode.cc,v
retrieving revision 1.36
diff -c -r1.36 lsode.cc
*** src/DLD-FUNCTIONS/lsode.cc 20 Mar 2006 18:29:20 -0000 1.36
--- src/DLD-FUNCTIONS/lsode.cc 2 Apr 2006 20:23:31 -0000
***************
*** 38,43 ****
--- 38,44 ----
#include "gripes.h"
#include "oct-obj.h"
#include "ov-fcn.h"
+ #include "ov-cell.h"
#include "pager.h"
#include "pr-output.h"
#include "unwind-prot.h"
***************
*** 191,199 ****
state of the system @var{x_0}, so that the first row of the output\n\
is @var{x_0}.\n\
\n\
! The first argument, @var{fcn}, is a string that names the function to\n\
! call to compute the vector of right hand sides for the set of equations.\n\
! The function must have the form\n\
\n\
@example\n\
@var{xdot} = f (@var{x}, @var{t})\n\
--- 192,201 ----
state of the system @var{x_0}, so that the first row of the output\n\
is @var{x_0}.\n\
\n\
! The first argument, @var{fcn}, is a string, or cell array of strings,\n\
! inline or function handles, that names the function to call to compute\n\
! the vector of right hand sides for the set of equations. The function\n\
! must have the form\n\
\n\
@example\n\
@var{xdot} = f (@var{x}, @var{t})\n\
***************
*** 286,333 ****
if (nargin > 2 && nargin < 5 && nargout < 4)
{
lsode_fcn = 0;
lsode_jac = 0;
octave_value f_arg = args(0);
! switch (f_arg.rows ())
{
! case 1:
! lsode_fcn = extract_function
! (f_arg, "lsode", "__lsode_fcn__",
! "function xdot = __lsode_fcn__ (x, t) xdot = ",
! "; endfunction");
! break;
!
! case 2:
! {
! string_vector tmp = f_arg.all_strings ();
!
! if (! error_state)
! {
! lsode_fcn = extract_function
! (tmp(0), "lsode", "__lsode_fcn__",
! "function xdot = __lsode_fcn__ (x, t) xdot = ",
! "; endfunction");
! if (lsode_fcn)
{
! lsode_jac = extract_function
! (tmp(1), "lsode", "__lsode_jac__",
! "function jac = __lsode_jac__ (x, t) jac = ",
! "; endfunction");
! if (! lsode_jac)
! lsode_fcn = 0;
}
! }
! }
! break;
!
! default:
! LSODE_ABORT1
! ("first arg should be a string or 2-element string array");
}
if (error_state || ! lsode_fcn)
--- 288,404 ----
if (nargin > 2 && nargin < 5 && nargout < 4)
{
+ std::string fcn_name, fname, jac_name, jname;
lsode_fcn = 0;
lsode_jac = 0;
octave_value f_arg = args(0);
! if (f_arg.is_cell ())
! {
! Cell c = f_arg.cell_value ();
! if (c.length() == 1)
! f_arg = c(0);
! else if (c.length() == 2)
! {
! if (c(0).is_function_handle () || c(0).is_inline_function ())
! lsode_fcn = c(0).function_value ();
! else
! {
! fcn_name = unique_symbol_name ("__lsode_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, t) y = ");
! lsode_fcn = extract_function
! (c(0), "lsode", fcn_name, fname, "; endfunction");
! }
!
! if (lsode_fcn)
! {
! if (c(1).is_function_handle () || c(1).is_inline_function ())
! lsode_jac = c(1).function_value ();
! else
! {
! jac_name = unique_symbol_name ("__lsode_jac__");
! jname = "function jac = ";
! jname.append(jac_name);
! jname.append (" (x, t) jac = ");
! lsode_jac = extract_function
! (c(1), "lsode", jac_name, jname, "; endfunction");
!
! if (!lsode_jac)
! {
! if (fcn_name.length())
! clear_function (fcn_name);
! lsode_fcn = 0;
! }
! }
! }
! }
! else
! LSODE_ABORT1 ("incorrect number of elements in cell array");
! }
!
! if (!lsode_fcn && ! f_arg.is_cell())
{
! if (f_arg.is_function_handle () || f_arg.is_inline_function ())
! lsode_fcn = f_arg.function_value ();
! else
! {
! switch (f_arg.rows ())
! {
! case 1:
! do
! {
! fcn_name = unique_symbol_name ("__lsode_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, t) y = ");
! lsode_fcn = extract_function
! (f_arg, "lsode", fcn_name, fname, "; endfunction");
! }
! while (0);
! break;
! case 2:
{
! string_vector tmp = f_arg.all_strings ();
! if (! error_state)
! {
! fcn_name = unique_symbol_name ("__lsode_fcn__");
! fname = "function y = ";
! fname.append (fcn_name);
! fname.append (" (x, t) y = ");
! lsode_fcn = extract_function
! (tmp(0), "lsode", fcn_name, fname, "; endfunction");
!
! if (lsode_fcn)
! {
! jac_name = unique_symbol_name ("__lsode_jac__");
! jname = "function jac = ";
! jname.append(jac_name);
! jname.append (" (x, t) jac = ");
! lsode_jac = extract_function
! (tmp(1), "lsode", jac_name, jname,
! "; endfunction");
!
! if (!lsode_jac)
! {
! if (fcn_name.length())
! clear_function (fcn_name);
! lsode_fcn = 0;
! }
! }
! }
}
! break;
!
! default:
! LSODE_ABORT1
! ("first arg should be a string or 2-element string array");
! }
! }
}
if (error_state || ! lsode_fcn)
***************
*** 372,377 ****
--- 443,453 ----
else
output = ode.integrate (out_times);
+ if (fcn_name.length())
+ clear_function (fcn_name);
+ if (jac_name.length())
+ clear_function (jac_name);
+
if (! error_state)
{
std::string msg = ode.error_message ();
Index: src/ls-hdf5.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ls-hdf5.h,v
retrieving revision 1.8
diff -c -r1.8 ls-hdf5.h
*** src/ls-hdf5.h 16 May 2005 20:07:36 -0000 1.8
--- src/ls-hdf5.h 1 Apr 2006 13:55:11 -0000
***************
*** 67,73 ****
}
}
! void open (const char *name, int mode, int prot = 0)
{
clear ();
--- 67,73 ----
}
}
! void open (const char *name, int mode, int)
{
clear ();
Index: test/test_number.m
===================================================================
RCS file: /usr/local/cvsroot/octave/test/test_number.m,v
retrieving revision 1.1
diff -c -r1.1 test_number.m
*** test/test_number.m 6 Jan 2006 00:24:06 -0000 1.1
--- test/test_number.m 2 Apr 2006 20:23:32 -0000
***************
*** 55,66 ****
%% test/octave.test/number/isvector-5.m
%!test
%! warn_str_to_num = 0;
! %! assert(!(isvector ("t")));
%% test/octave.test/number/isvector-6.m
%!test
%! warn_str_to_num = 0;
! %! assert(!(isvector ("test")));
%% test/octave.test/number/isvector-7.m
%!assert(!(isvector (["test"; "ing"])));
--- 55,66 ----
%% test/octave.test/number/isvector-5.m
%!test
%! warn_str_to_num = 0;
! %! assert((isvector ("t")));
%% test/octave.test/number/isvector-6.m
%!test
%! warn_str_to_num = 0;
! %! assert((isvector ("test")));
%% test/octave.test/number/isvector-7.m
%!assert(!(isvector (["test"; "ing"])));
***************
*** 68,74 ****
%% test/octave.test/number/isvector-8.m
%!test
%! s.a = 1;
! %! assert(!(isvector (s)));
%% test/octave.test/number/isvector-9.m
%!error isvector ();
--- 68,74 ----
%% test/octave.test/number/isvector-8.m
%!test
%! s.a = 1;
! %! assert((isvector (s)));
%% test/octave.test/number/isvector-9.m
%!error isvector ();
***************
*** 91,97 ****
%% test/octave.test/number/isscalar-5.m
%!test
%! warn_str_to_num = 0;
! %! assert(!(isscalar ("t")));
%% test/octave.test/number/isscalar-6.m
%!assert(!(isscalar ("test")));
--- 91,97 ----
%% test/octave.test/number/isscalar-5.m
%!test
%! warn_str_to_num = 0;
! %! assert((isscalar ("t")));
%% test/octave.test/number/isscalar-6.m
%!assert(!(isscalar ("test")));
***************
*** 102,108 ****
%% test/octave.test/number/isscalar-8.m
%!test
%! s.a = 1;
! %! assert(!(isscalar (s)));
%% test/octave.test/number/isscalar-9.m
%!error isscalar ();
--- 102,108 ----
%% test/octave.test/number/isscalar-8.m
%!test
%! s.a = 1;
! %! assert((isscalar (s)));
%% test/octave.test/number/isscalar-9.m
%!error isscalar ();
Index: doc/interpreter/matrix.txi
===================================================================
RCS file: /usr/local/cvsroot/octave/doc/interpreter/matrix.txi,v
retrieving revision 1.12
diff -c -r1.12 matrix.txi
*** doc/interpreter/matrix.txi 3 Jun 2004 19:32:02 -0000 1.12
--- doc/interpreter/matrix.txi 2 Apr 2006 20:23:24 -0000
***************
*** 143,150 ****
@DOCSTRING(randn)
! The @code{rand} and @code{randn} functions use separate generators.
! This ensures that
@example
@group
--- 143,157 ----
@DOCSTRING(randn)
! @DOCSTRING(rande)
!
! @DOCSTRING(randp)
!
! @DOCSTRING(randg)
!
! The new random generators all use a common Mersenne Twister generator,
! and so the state of only one of the generators needs to be reset.
! The old generator function use separate generators. This ensures that
@example
@group
Index: libcruft/ranlib/wrap.f
===================================================================
RCS file: /usr/local/cvsroot/octave/libcruft/ranlib/wrap.f,v
retrieving revision 1.1
diff -c -r1.1 wrap.f
*** libcruft/ranlib/wrap.f 19 Jul 1996 01:29:53 -0000 1.1
--- libcruft/ranlib/wrap.f 2 Apr 2006 20:23:25 -0000
***************
*** 8,10 ****
--- 8,25 ----
result = genunf (real (low), real (high))
return
end
+ subroutine dgenexp (av, result)
+ double precision av, result
+ result = genexp (real (av))
+ return
+ end
+ subroutine dgengam (a, r, result)
+ double precision a, r, result
+ result = gengam (real (a), real (r))
+ return
+ end
+ subroutine dignpoi (mu, result)
+ double precision mu, result
+ result = ignpoi (real (mu))
+ return
+ end
Index: liboctave/Makefile.in
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/Makefile.in,v
retrieving revision 1.213
diff -c -r1.213 Makefile.in
*** liboctave/Makefile.in 16 Mar 2006 17:48:56 -0000 1.213
--- liboctave/Makefile.in 2 Apr 2006 20:23:25 -0000
***************
*** 69,76 ****
oct-passwd.h oct-rand.h oct-rl-edit.h oct-rl-hist.h \
oct-shlib.h oct-sort.h oct-spparms.h oct-syscalls.h \
oct-sparse.h oct-time.h oct-types.h oct-uname.h \
! pathlen.h pathsearch.h \
! prog-args.h so-array.h sparse-sort.h statdefs.h str-vec.h \
sparse-util.h sun-utils.h sysdir.h systime.h syswait.h \
$(OPTS_INC) \
$(MATRIX_INC) \
--- 69,77 ----
oct-passwd.h oct-rand.h oct-rl-edit.h oct-rl-hist.h \
oct-shlib.h oct-sort.h oct-spparms.h oct-syscalls.h \
oct-sparse.h oct-time.h oct-types.h oct-uname.h \
! pathlen.h pathsearch.h prog-args.h \
! randgamma.h randmtzig.h randpoisson.h \
! so-array.h sparse-sort.h statdefs.h str-vec.h \
sparse-util.h sun-utils.h sysdir.h systime.h syswait.h \
$(OPTS_INC) \
$(MATRIX_INC) \
***************
*** 128,134 ****
$(SPARSE_MX_OP_SRC)
LIBOCTAVE_C_SOURCES := f2c-main.c filemode.c getopt.c getopt1.c \
! lo-cieee.c lo-cutils.c mkdir.c oct-getopt.c rename.c \
rmdir.c strftime.c strptime.c tempname.c tempnam.c
LIBOCTAVE_SOURCES := $(LIBOCTAVE_CXX_SOURCES) $(LIBOCTAVE_C_SOURCES)
--- 129,136 ----
$(SPARSE_MX_OP_SRC)
LIBOCTAVE_C_SOURCES := f2c-main.c filemode.c getopt.c getopt1.c \
! lo-cieee.c lo-cutils.c mkdir.c oct-getopt.c \
! randgamma.c randmtzig.c randpoisson.c rename.c \
rmdir.c strftime.c strptime.c tempname.c tempnam.c
LIBOCTAVE_SOURCES := $(LIBOCTAVE_CXX_SOURCES) $(LIBOCTAVE_C_SOURCES)
Index: liboctave/oct-rand.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/oct-rand.cc,v
retrieving revision 1.6
diff -c -r1.6 oct-rand.cc
*** liboctave/oct-rand.cc 26 Apr 2005 19:24:29 -0000 1.6
--- liboctave/oct-rand.cc 2 Apr 2006 20:23:26 -0000
***************
*** 24,45 ****
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include "f77-fcn.h"
#include "lo-error.h"
#include "oct-rand.h"
#include "oct-time.h"
// Possible distributions of random numbers. This was handled with an
// enum, but unwind_protecting that doesn't work so well.
#define uniform_dist 1
#define normal_dist 2
// Current distribution of random numbers.
static int current_distribution = uniform_dist;
! // Has the seed been set yet?
! static bool initialized = false;
extern "C"
{
--- 24,57 ----
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
+ #include <vector>
#include "f77-fcn.h"
+ #include "lo-ieee.h"
#include "lo-error.h"
+ #include "lo-mappers.h"
#include "oct-rand.h"
#include "oct-time.h"
+ #include "data-conv.h"
+ #include "randmtzig.h"
+ #include "randpoisson.h"
+ #include "randgamma.h"
// Possible distributions of random numbers. This was handled with an
// enum, but unwind_protecting that doesn't work so well.
#define uniform_dist 1
#define normal_dist 2
+ #define expon_dist 3
+ #define poisson_dist 4
+ #define gamma_dist 5
// Current distribution of random numbers.
static int current_distribution = uniform_dist;
! // Has the seed/state been set yet?
! static bool old_initialized = false;
! static bool new_initialized = false;
! static bool use_old_generators = false;
extern "C"
{
***************
*** 50,55 ****
--- 62,76 ----
F77_FUNC (dgenunf, DGENUNF) (const double&, const double&, double&);
F77_RET_T
+ F77_FUNC (dgenexp, DGENEXP) (const double&, double&);
+
+ F77_RET_T
+ F77_FUNC (dignpoi, DIGNPOI) (const double&, double&);
+
+ F77_RET_T
+ F77_FUNC (dgengam, DGENGAM) (const double&, const double&, double&);
+
+ F77_RET_T
F77_FUNC (setall, SETALL) (const octave_idx_type&, const octave_idx_type&);
F77_RET_T
***************
*** 83,89 ****
// work ok to give fairly different seeds each time Octave starts.
static void
! do_initialization (void)
{
octave_localtime tm;
--- 104,110 ----
// work ok to give fairly different seeds each time Octave starts.
static void
! do_old_initialization (void)
{
octave_localtime tm;
***************
*** 99,118 ****
F77_FUNC (setall, SETALL) (s0, s1);
! initialized = true;
}
static inline void
maybe_initialize (void)
{
! if (! initialized)
! do_initialization ();
}
double
octave_rand::seed (void)
{
! maybe_initialize ();
union d2i { double d; octave_idx_type i[2]; };
union d2i u;
--- 120,151 ----
F77_FUNC (setall, SETALL) (s0, s1);
! old_initialized = true;
}
static inline void
maybe_initialize (void)
{
! if (use_old_generators)
! {
! if (! old_initialized)
! do_old_initialization ();
! }
! else
! {
! if (! new_initialized)
! {
! oct_init_by_entropy ();
! new_initialized = true;
! }
! }
}
double
octave_rand::seed (void)
{
! if (! old_initialized)
! do_old_initialization ();
union d2i { double d; octave_idx_type i[2]; };
union d2i u;
***************
*** 123,128 ****
--- 156,162 ----
void
octave_rand::seed (double s)
{
+ use_old_generators = true;
maybe_initialize ();
union d2i { double d; octave_idx_type i[2]; };
***************
*** 133,138 ****
--- 167,207 ----
F77_FUNC (setsd, SETSD) (i0, i1);
}
+ ColumnVector
+ octave_rand::state (void)
+ {
+ ColumnVector s (MT_N + 1);
+ if (! new_initialized)
+ {
+ oct_init_by_entropy ();
+ new_initialized = true;
+ }
+
+ OCTAVE_LOCAL_BUFFER (unsigned FOUR_BYTE_INT, tmp, MT_N + 1);
+ oct_get_state (tmp);
+ for (octave_idx_type i = 0; i <= MT_N; i++)
+ s.elem (i) = static_cast<double>(tmp [i]);
+ return s;
+ }
+
+ void
+ octave_rand::state (const ColumnVector &s)
+ {
+ use_old_generators = false;
+ maybe_initialize ();
+
+ octave_idx_type len = s.length();
+ octave_idx_type n = len < MT_N + 1 ? len : MT_N + 1;
+ OCTAVE_LOCAL_BUFFER (unsigned FOUR_BYTE_INT, tmp, MT_N + 1);
+ for (octave_idx_type i = 0; i < n; i++)
+ tmp[i] = static_cast<unsigned FOUR_BYTE_INT> (s.elem(i));
+
+ if (len == MT_N + 1 && tmp[MT_N] <= MT_N && tmp[MT_N] > 0)
+ oct_set_state (tmp);
+ else
+ oct_init_by_array (tmp, len);
+ }
+
std::string
octave_rand::distribution (void)
{
***************
*** 142,147 ****
--- 211,222 ----
return "uniform";
else if (current_distribution == normal_dist)
return "normal";
+ else if (current_distribution == expon_dist)
+ return "exponential";
+ else if (current_distribution == poisson_dist)
+ return "poisson";
+ else if (current_distribution == gamma_dist)
+ return "gamma";
else
{
abort ();
***************
*** 156,161 ****
--- 231,242 ----
octave_rand::uniform_distribution ();
else if (d == "normal")
octave_rand::normal_distribution ();
+ else if (d == "exponential")
+ octave_rand::exponential_distribution ();
+ else if (d == "poisson")
+ octave_rand::poisson_distribution ();
+ else if (d == "gamma")
+ octave_rand::gamma_distribution ();
else
(*current_liboctave_error_handler) ("rand: invalid distribution");
}
***************
*** 180,319 ****
F77_FUNC (setcgn, SETCGN) (normal_dist);
}
! double
! octave_rand::scalar (void)
{
maybe_initialize ();
! double retval = 0.0;
! switch (current_distribution)
! {
! case uniform_dist:
! F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval);
! break;
! case normal_dist:
! F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval);
! break;
! default:
! abort ();
! break;
! }
! return retval;
}
! #define MAKE_RAND_MAT(mat, nr, nc, f, F) \
! do \
! { \
! double val; \
! for (volatile octave_idx_type j = 0; j < nc; j++) \
! for (volatile octave_idx_type i = 0; i < nr; i++) \
! { \
! OCTAVE_QUIT; \
! F77_FUNC (f, F) (0.0, 1.0, val); \
! mat(i,j) = val; \
! } \
! } \
! while (0)
! Matrix
! octave_rand::matrix (octave_idx_type n, octave_idx_type m)
{
maybe_initialize ();
! Matrix retval;
! if (n >= 0 && m >= 0)
{
! retval.resize (n, m);
!
! if (n > 0 && m > 0)
{
! switch (current_distribution)
{
! case uniform_dist:
! MAKE_RAND_MAT (retval, n, m, dgenunf, DGENUNF);
! break;
!
! case normal_dist:
! MAKE_RAND_MAT (retval, n, m, dgennor, DGENNOR);
! break;
!
! default:
! abort ();
! break;
}
}
}
else
! (*current_liboctave_error_handler) ("rand: invalid negative argument");
return retval;
}
! #define MAKE_RAND_ND_ARRAY(mat, len, f, F) \
do \
{ \
double val; \
for (volatile octave_idx_type i = 0; i < len; i++) \
{ \
OCTAVE_QUIT; \
! F77_FUNC (f, F) (0.0, 1.0, val); \
! mat(i) = val; \
} \
} \
while (0)
! NDArray
! octave_rand::nd_array (const dim_vector& dims)
{
maybe_initialize ();
! NDArray retval;
! if (! dims.all_zero ())
{
! retval.resize (dims);
! octave_idx_type len = retval.length ();
! switch (current_distribution)
{
! case uniform_dist:
! MAKE_RAND_ND_ARRAY (retval, len, dgenunf, DGENUNF);
! break;
! case normal_dist:
! MAKE_RAND_ND_ARRAY (retval, len, dgennor, DGENNOR);
! break;
! default:
! abort ();
! break;
}
}
return retval;
}
! #define MAKE_RAND_ARRAY(array, n, f, F) \
! do \
! { \
! double val; \
! for (volatile octave_idx_type i = 0; i < n; i++) \
! { \
! OCTAVE_QUIT; \
! F77_FUNC (f, F) (0.0, 1.0, val); \
! array(i) = val; \
! } \
! } \
! while (0)
Array<double>
! octave_rand::vector (octave_idx_type n)
{
maybe_initialize ();
--- 261,512 ----
F77_FUNC (setcgn, SETCGN) (normal_dist);
}
! void
! octave_rand::exponential_distribution (void)
{
maybe_initialize ();
! current_distribution = expon_dist;
! F77_FUNC (setcgn, SETCGN) (expon_dist);
! }
! void
! octave_rand::poisson_distribution (void)
! {
! maybe_initialize ();
! current_distribution = poisson_dist;
! F77_FUNC (setcgn, SETCGN) (poisson_dist);
}
! void
! octave_rand::gamma_distribution (void)
! {
! maybe_initialize ();
! current_distribution = gamma_dist;
!
! F77_FUNC (setcgn, SETCGN) (gamma_dist);
! }
!
!
! double
! octave_rand::scalar (double a)
{
maybe_initialize ();
! double retval = 0.0;
! if (use_old_generators)
{
! switch (current_distribution)
{
! case uniform_dist:
! F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval);
! break;
!
! case normal_dist:
! F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval);
! break;
!
! case expon_dist:
! F77_FUNC (dgenexp, DGENEXP) (1.0, retval);
! break;
!
! case poisson_dist:
! if (a < 0.0 || xisnan(a) || xisinf(a))
! retval = octave_NaN;
! else
{
! // workaround bug in ignpoi, by calling with different Mu
! F77_FUNC (dignpoi, DIGNPOI) (a + 1, retval);
! F77_FUNC (dignpoi, DIGNPOI) (a, retval);
}
+ break;
+
+ case gamma_dist:
+ if (a <= 0.0 || xisnan(a) || xisinf(a))
+ retval = octave_NaN;
+ else
+ F77_FUNC (dgengam, DGENGAM) (1.0, a, retval);
+ break;
+
+ default:
+ abort ();
+ break;
}
}
else
! {
! switch (current_distribution)
! {
! case uniform_dist:
! retval = oct_randu();
! break;
!
! case normal_dist:
! retval = oct_randn();
! break;
!
! case expon_dist:
! retval = oct_rande();
! break;
!
! case poisson_dist:
! retval = oct_randp(a);
! break;
!
! case gamma_dist:
! retval = oct_randg(a);
! break;
!
! default:
! abort ();
! break;
! }
! }
return retval;
}
! #define MAKE_RAND(len) \
do \
{ \
double val; \
for (volatile octave_idx_type i = 0; i < len; i++) \
{ \
OCTAVE_QUIT; \
! RAND_FUNC (val); \
! v[i] = val; \
} \
} \
while (0)
! static void
! fill_rand (octave_idx_type len, double *v, double a)
{
maybe_initialize ();
! if (len < 1)
! return;
! switch (current_distribution)
{
! case uniform_dist:
! if (use_old_generators)
! {
! #define RAND_FUNC(x) F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, x)
! MAKE_RAND (len);
! #undef RAND_FUNC
! }
! else
! oct_fill_randu (len, v);
! break;
! case normal_dist:
! if (use_old_generators)
! {
! #define RAND_FUNC(x) F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, x)
! MAKE_RAND (len);
! #undef RAND_FUNC
! }
! else
! oct_fill_randn (len, v);
! break;
! case expon_dist:
! if (use_old_generators)
{
! #define RAND_FUNC(x) F77_FUNC (dgenexp, DGENEXP) (1.0, x)
! MAKE_RAND (len);
! #undef RAND_FUNC
! }
! else
! oct_fill_rande (len, v);
! break;
! case poisson_dist:
! if (use_old_generators)
! {
! if (a < 0.0 || xisnan(a) || xisinf(a))
! #define RAND_FUNC(x) x = octave_NaN;
! MAKE_RAND (len);
! #undef RAND_FUNC
! else
! {
! // workaround bug in ignpoi, by calling with different Mu
! double tmp;
! F77_FUNC (dignpoi, DIGNPOI) (a + 1, tmp);
! #define RAND_FUNC(x) F77_FUNC (dignpoi, DIGNPOI) (a, x)
! MAKE_RAND (len);
! #undef RAND_FUNC
! }
! }
! else
! oct_fill_randp (a, len, v);
! break;
! case gamma_dist:
! if (use_old_generators)
! {
! if (a <= 0.0 || xisnan(a) || xisinf(a))
! #define RAND_FUNC(x) x = octave_NaN;
! MAKE_RAND (len);
! #undef RAND_FUNC
! else
! #define RAND_FUNC(x) F77_FUNC (dgengam, DGENGAM) (1.0, a, x)
! MAKE_RAND (len);
! #undef RAND_FUNC
}
+ else
+ oct_fill_randg (a, len, v);
+ break;
+
+ default:
+ abort ();
+ break;
+ }
+
+ return;
+ }
+
+ Matrix
+ octave_rand::matrix (octave_idx_type n, octave_idx_type m, double a)
+ {
+ Matrix retval;
+
+ if (n >= 0 && m >= 0)
+ {
+ retval.resize (n, m);
+
+ if (n > 0 && m > 0)
+ fill_rand (retval.capacity(), retval.fortran_vec(), a);
}
+ else
+ (*current_liboctave_error_handler) ("rand: invalid negative argument");
return retval;
}
! NDArray
! octave_rand::nd_array (const dim_vector& dims, double a)
! {
! NDArray retval;
!
! if (! dims.all_zero ())
! {
! retval.resize (dims);
!
! fill_rand (retval.capacity(), retval.fortran_vec(), a);
! }
!
! return retval;
! }
Array<double>
! octave_rand::vector (octave_idx_type n, double a)
{
maybe_initialize ();
***************
*** 323,342 ****
{
retval.resize (n);
! switch (current_distribution)
! {
! case uniform_dist:
! MAKE_RAND_ARRAY (retval, n, dgenunf, DGENUNF);
! break;
!
! case normal_dist:
! MAKE_RAND_ARRAY (retval, n, dgennor, DGENNOR);
! break;
!
! default:
! abort ();
! break;
! }
}
else if (n < 0)
(*current_liboctave_error_handler) ("rand: invalid negative argument");
--- 516,522 ----
{
retval.resize (n);
! fill_rand (retval.capacity(), retval.fortran_vec(), a);
}
else if (n < 0)
(*current_liboctave_error_handler) ("rand: invalid negative argument");
Index: liboctave/oct-rand.h
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/oct-rand.h,v
retrieving revision 1.4
diff -c -r1.4 oct-rand.h
*** liboctave/oct-rand.h 26 Apr 2005 19:24:29 -0000 1.4
--- liboctave/oct-rand.h 2 Apr 2006 20:23:26 -0000
***************
*** 26,31 ****
--- 26,32 ----
#include <string>
+ #include "dColVector.h"
#include "dMatrix.h"
#include "dNDArray.h"
***************
*** 38,43 ****
--- 39,50 ----
// Set the seed.
static void seed (double s);
+ // Return the current state.
+ static ColumnVector state (void);
+
+ // Set the current state/
+ static void state (const ColumnVector &s);
+
// Return the current distribution.
static std::string distribution (void);
***************
*** 49,67 ****
static void normal_distribution (void);
// Return the next number from the sequence.
! static double scalar (void);
// Return a matrix of numbers from the sequence, filled in column
// major order.
! static Matrix matrix (octave_idx_type r, octave_idx_type c);
// Return an N-dimensional array of numbers from the sequence,
// filled in column major order.
! static NDArray nd_array (const dim_vector& dims);
// Return an array of numbers from the sequence.
! static Array<double> vector (octave_idx_type n);
};
#endif
--- 56,80 ----
static void normal_distribution (void);
+ static void exponential_distribution (void);
+
+ static void poisson_distribution (void);
+
+ static void gamma_distribution (void);
+
// Return the next number from the sequence.
! static double scalar (double a = 1.);
// Return a matrix of numbers from the sequence, filled in column
// major order.
! static Matrix matrix (octave_idx_type r, octave_idx_type c, double a = 1.);
// Return an N-dimensional array of numbers from the sequence,
// filled in column major order.
! static NDArray nd_array (const dim_vector& dims, double a = 1.);
// Return an array of numbers from the sequence.
! static Array<double> vector (octave_idx_type n, double a = 1.);
};
#endif
Index: src/DLD-FUNCTIONS/rand.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/DLD-FUNCTIONS/rand.cc,v
retrieving revision 1.26
diff -c -r1.26 rand.cc
*** src/DLD-FUNCTIONS/rand.cc 18 May 2005 11:43:38 -0000 1.26
--- src/DLD-FUNCTIONS/rand.cc 2 Apr 2006 20:23:31 -0000
***************
*** 42,69 ****
#include "utils.h"
static octave_value
! do_rand (const octave_value_list& args, int nargin, const char *fcn)
{
octave_value retval;
!
dim_vector dims;
switch (nargin)
{
case 0:
{
! dims.resize (2);
!
! dims(0) = 1;
! dims(1) = 1;
goto gen_matrix;
}
break;
case 1:
{
! octave_value tmp = args(0);
if (tmp.is_string ())
{
--- 42,97 ----
#include "utils.h"
static octave_value
! do_rand (const octave_value_list& args, int nargin, const char *fcn,
! bool additional_arg = false)
{
octave_value retval;
! NDArray a;
! int idx = 0;
dim_vector dims;
+ if (additional_arg)
+ {
+ if (nargin == 0)
+ {
+ error ("%s: expecting at least one argument", fcn);
+ goto done;
+ }
+ else if (args(0).is_string())
+ additional_arg = false;
+ else
+ {
+ a = args(0).array_value ();
+ if (error_state)
+ {
+ error ("%s: expecting scalar or matrix arguments", fcn);
+ goto done;
+ }
+ idx++;
+ nargin--;
+ }
+ }
+
switch (nargin)
{
case 0:
{
! if (additional_arg)
! dims = a.dims ();
! else
! {
! dims.resize (2);
+ dims(0) = 1;
+ dims(1) = 1;
+ }
goto gen_matrix;
}
break;
case 1:
{
! octave_value tmp = args(idx);
if (tmp.is_string ())
{
***************
*** 73,82 ****
{
retval = octave_rand::distribution ();
}
! else if (s_arg == "seed" || s_arg == "state")
{
retval = octave_rand::seed ();
}
else if (s_arg == "uniform")
{
octave_rand::uniform_distribution ();
--- 101,114 ----
{
retval = octave_rand::distribution ();
}
! else if (s_arg == "seed")
{
retval = octave_rand::seed ();
}
+ else if (s_arg == "state" || s_arg == "twister")
+ {
+ retval = octave_rand::state ();
+ }
else if (s_arg == "uniform")
{
octave_rand::uniform_distribution ();
***************
*** 85,90 ****
--- 117,134 ----
{
octave_rand::normal_distribution ();
}
+ else if (s_arg == "exponential")
+ {
+ octave_rand::exponential_distribution ();
+ }
+ else if (s_arg == "poisson")
+ {
+ octave_rand::poisson_distribution ();
+ }
+ else if (s_arg == "gamma")
+ {
+ octave_rand::gamma_distribution ();
+ }
else
error ("%s: unrecognized string argument", fcn);
}
***************
*** 176,194 ****
default:
{
! octave_value tmp = args(0);
if (nargin == 2 && tmp.is_string ())
{
std::string ts = tmp.string_value ();
! if (ts == "seed" || ts == "state")
{
! double d = args(1).double_value ();
if (! error_state)
octave_rand::seed (d);
}
else
error ("%s: unrecognized string argument", fcn);
}
--- 220,246 ----
default:
{
! octave_value tmp = args(idx);
if (nargin == 2 && tmp.is_string ())
{
std::string ts = tmp.string_value ();
! if (ts == "seed")
{
! double d = args(idx+1).double_value ();
if (! error_state)
octave_rand::seed (d);
}
+ else if (ts == "state" || ts == "twister")
+ {
+ ColumnVector s =
+ ColumnVector (args(idx+1).vector_value(false, true));
+
+ if (! error_state)
+ octave_rand::state (s);
+ }
else
error ("%s: unrecognized string argument", fcn);
}
***************
*** 198,204 ****
for (int i = 0; i < nargin; i++)
{
! dims(i) = (octave_idx_type)args(i).int_value ();
if (error_state)
{
--- 250,256 ----
for (int i = 0; i < nargin; i++)
{
! dims(i) = (octave_idx_type)args(idx+i).int_value ();
if (error_state)
{
***************
*** 221,253 ****
dims.chop_trailing_singletons ();
! return octave_rand::nd_array (dims);
}
DEFUN_DLD (rand, args, ,
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} rand (@var{x})\n\
@deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\
! @deftypefnx {Loadable Function} {} rand (@code{\"seed\"}, @var{x})\n\
Return a matrix with random elements uniformly distributed on the\n\
interval (0, 1). The arguments are handled the same as the arguments\n\
! for @code{eye}. In\n\
! addition, you can set the seed for the random number generator using the\n\
form\n\
\n\
@example\n\
! rand (\"seed\", @var{x})\n\
@end example\n\
\n\
! @noindent\n\
! where @var{x} is a scalar value. If called as\n\
\n\
@example\n\
! rand (\"seed\")\n\
@end example\n\
\n\
@noindent\n\
! @code{rand} returns the current value of the seed.\n\
@end deftypefn")
{
octave_value retval;
--- 273,369 ----
dims.chop_trailing_singletons ();
! if (additional_arg)
! {
! if (a.length() == 1)
! return octave_rand::nd_array (dims, a(0));
! else
! {
! if (a.dims() != dims)
! {
! error ("%s: mismatch in argument size", fcn);
! return retval;
! }
! octave_idx_type len = a.length ();
! NDArray m (dims);
! double *v = m.fortran_vec ();
! for (octave_idx_type i = 0; i < len; i++)
! v[i] = octave_rand::scalar (a(i));
! return m;
! }
! }
! else
! return octave_rand::nd_array (dims);
}
DEFUN_DLD (rand, args, ,
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} rand (@var{x})\n\
@deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\
! @deftypefnx {Loadable Function} {} rand (\"state\", @var{x})\n\
! @deftypefnx {Loadable Function} {} rand (\"seed\", @var{x})\n\
Return a matrix with random elements uniformly distributed on the\n\
interval (0, 1). The arguments are handled the same as the arguments\n\
! for @code{eye}.\n\
! \n\
! You can query the state of the random number generator using the\n\
form\n\
\n\
@example\n\
! v = rand (\"state\")\n\
@end example\n\
\n\
! This returns a column vector @var{v} of length 625. Later, you can\n\
! restore the random number generator to the state @var{v}\n\
! using the form\n\
\n\
@example\n\
! rand (\"state\", v)\n\
@end example\n\
\n\
@noindent\n\
! You may also initialize the state vector from an arbitrary vector of\n\
! length <= 625 for @var{v}. This new state will be a hash based on the\n\
! the value of @var{v}, not @var{v} itself.\n\
! \n\
! By default, the generator is initialized from @code{/dev/urandom} if it is\n\
! available, otherwise from cpu time, wall clock time and the current\n\
! fraction of a second.\n\
! \n\
! @code{rand} uses the Mersenne Twister with a period of 2^19937-1\n\
! (See M. Matsumoto and T. Nishimura, ``Mersenne Twister: A 623-dimensionally\n\
! equidistributed uniform pseudorandom number generator'', ACM Trans. on\n\
! Modeling and Computer Simulation Vol. 8, No. 1, Januray pp.3-30 1998,\n\
! @url{http://www.math.keio.ac.jp/~matumoto/emt.html}).\n\
! Do NOT use for CRYPTOGRAPHY without securely hashing several returned\n\
! values together, otherwise the generator state can be learned after\n\
! reading 624 consecutive values.\n\
! \n\
! @code{rand} includes a second random number generator, that was the\n\
! previous generator used in octave. The new generator is used by default\n\
! as it is significantly faster than the old generator, and produces\n\
! random numebrs with a significantly longer cycle time. However, in\n\
! some circumstances it might be desireable to obtain the same random\n\
! sequences as used by the old generators. To do this the keyword\n\
! \"seed\" is used to specify that the old generators should be use,\n\
! as in\n\
! \n\
! @example\n\
! rand (\"seed\", val)\n\
! @end example\n\
! \n\
! which sets the seed of the generator to @var{val}. The seed of the\n\
! generator can be queried with\n\
! \n\
! @example\n\
! s = rand (\"seed\")\n\
! @end example\n\
! \n\
! However, it should be noted that querying the seed will not cause\n\
! @code{rand} to use the old generators, only setting the seed will.\n\
! To cause @code{rand} to once again use the new generators, the\n\
! keyword \"state\" should be used to reset the state of the @code{rand}.\n\
! @seealso{randn,rande,randg,randp}\n\
@end deftypefn")
{
octave_value retval;
***************
*** 259,264 ****
--- 375,438 ----
return retval;
}
+ /*
+ %!test # 'state' can be a scalar
+ %! rand('state',12); x = rand(1,4);
+ %! rand('state',12); y = rand(1,4);
+ %! assert(x,y);
+ %!test # 'state' can be a vector
+ %! rand('state',[12,13]); x=rand(1,4);
+ %! rand('state',[12;13]); y=rand(1,4);
+ %! assert(x,y);
+ %!test # querying 'state' doesn't disturb sequence
+ %! rand('state',12); rand(1,2); x=rand(1,2);
+ %! rand('state',12); rand(1,2);
+ %! s=rand('state'); y=rand(1,2);
+ %! assert(x,y);
+ %! rand('state',s); z=rand(1,2);
+ %! assert(x,z);
+ %!test # 'seed' must be a scalar
+ %! rand('seed',12); x = rand(1,4);
+ %! rand('seed',12); y = rand(1,4);
+ %! assert(x,y);
+ %!error(rand('seed',[12,13]))
+ %!test # querying 'seed' returns a value which can be used later
+ %! s=rand('seed'); x=rand(1,2);
+ %! rand('seed',s); y=rand(1,2);
+ %! assert(x,y);
+ %!test # querying 'seed' doesn't disturb sequence
+ %! rand('seed',12); rand(1,2); x=rand(1,2);
+ %! rand('seed',12); rand(1,2);
+ %! s=rand('seed'); y=rand(1,2);
+ %! assert(x,y);
+ %! rand('seed',s); z=rand(1,2);
+ %! assert(x,z);
+ */
+
+ /*
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! rand("state",12);
+ %! x = rand(100000,1);
+ %! assert(max(x)<1.); %*** Please report this!!! ***
+ %! assert(min(x)>0.); %*** Please report this!!! ***
+ %! assert(mean(x),0.5,0.0024);
+ %! assert(var(x),1/48,0.0632);
+ %! assert(skewness(x),0,0.012);
+ %! assert(kurtosis(x),-6/5,0.0094);
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! rand("seed",12);
+ %! x = rand(100000,1);
+ %! assert(max(x)<1.); %*** Please report this!!! ***
+ %! assert(min(x)>0.); %*** Please report this!!! ***
+ %! assert(mean(x),0.5,0.0024);
+ %! assert(var(x),1/48,0.0632);
+ %! assert(skewness(x),0,0.012);
+ %! assert(kurtosis(x),-6/5,0.0094);
+ */
+
+
static std::string current_distribution = octave_rand::distribution ();
static void
***************
*** 271,295 ****
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} randn (@var{x})\n\
@deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\
! @deftypefnx {Loadable Function} {} randn (@code{\"seed\"}, @var{x})\n\
! Return a matrix with normally distributed random elements. The\n\
! arguments are handled the same as the arguments for @code{eye}. In\n\
! addition, you can set the seed for the random number generator using the\n\
! form\n\
\n\
! @example\n\
! randn (\"seed\", @var{x})\n\
! @end example\n\
\n\
! @noindent\n\
! where @var{x} is a scalar value. If called as\n\
! \n\
! @example\n\
! randn (\"seed\")\n\
! @end example\n\
! \n\
! @noindent\n\
! @code{randn} returns the current value of the seed.\n\
@end deftypefn")
{
octave_value retval;
--- 445,462 ----
"-*- texinfo -*-\n\
@deftypefn {Loadable Function} {} randn (@var{x})\n\
@deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\
! @deftypefnx {Loadable Function} {} randn (\"state\", @var{x})\n\
! @deftypefnx {Loadable Function} {} randn (\"seed\", @var{x})\n\
! Return a matrix with normally distributed random elements. The\n\
! arguments are handled the same as the arguments for @code{rand}.\n\
\n\
! By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\
! transform from a uniform to a normal distribution. (G. Marsaglia and\n\
! W.W. Tsang, 'Ziggurat method for generating random variables',\n\
! J. Statistical Software, vol 5, 2000,\n\
! @url{http://www.jstatsoft.org/v05/i08/})\n\
\n\
! @seealso{rand,rande,randg,randp}\n\
@end deftypefn")
{
octave_value retval;
***************
*** 318,323 ****
--- 485,853 ----
}
/*
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! rand("state",12);
+ %! x = randn(100000,1);
+ %! assert(mean(x),0,0.01);
+ %! assert(var(x),1,0.02);
+ %! assert(skewness(x),0,0.02);
+ %! assert(kurtosis(x),0,0.04);
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! rand("seed",12);
+ %! x = randn(100000,1);
+ %! assert(mean(x),0,0.01);
+ %! assert(var(x),1,0.02);
+ %! assert(skewness(x),0,0.02);
+ %! assert(kurtosis(x),0,0.04);
+ */
+
+ DEFUN_DLD (rande, args, ,
+ "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} rande (@var{x})\n\
+ @deftypefnx {Loadable Function} {} rande (@var{n}, @var{m})\n\
+ @deftypefnx {Loadable Function} {} rande (\"state\", @var{x})\n\
+ @deftypefnx {Loadable Function} {} rande (\"seed\", @var{x})\n\
+ Return a matrix with exponentially distributed random elements. The\n\
+ arguments are handled the same as the arguments for @code{rand}.\n\
+ \n\
+ By default, @code{randn} uses a Marsaglia and Tsang Ziggurat technique to\n\
+ transform from a uniform to a exponential distribution. (G. Marsaglia and\n\
+ W.W. Tsang, 'Ziggurat method for generating random variables',\n\
+ J. Statistical Software, vol 5, 2000,\n\
+ @url{http://www.jstatsoft.org/v05/i08/})\n\
+ @seealso{rand,randn,randg,randp}\n\
+ @end deftypefn")
+ {
+ octave_value retval;
+
+ int nargin = args.length ();
+
+ unwind_protect::begin_frame ("rande");
+
+ // This relies on the fact that elements are popped from the unwind
+ // stack in the reverse of the order they are pushed
+ // (i.e. current_distribution will be reset before calling
+ // reset_rand_generator()).
+
+ unwind_protect::add (reset_rand_generator, 0);
+ unwind_protect_str (current_distribution);
+
+ current_distribution = "exponential";
+
+ octave_rand::distribution (current_distribution);
+
+ retval = do_rand (args, nargin, "rande");
+
+ unwind_protect::run_frame ("rande");
+
+ return retval;
+ }
+
+ /*
+ %!test
+ %! % statistical tests may fail occasionally
+ %! rand("state",12);
+ %! x = rande(100000,1);
+ %! assert(min(x)>0); % *** Please report this!!! ***
+ %! assert(mean(x),1,0.01);
+ %! assert(var(x),1,0.03);
+ %! assert(skewness(x),2,0.06);
+ %! assert(kurtosis(x),6,0.7);
+ %!test
+ %! % statistical tests may fail occasionally
+ %! rand("seed",12);
+ %! x = rande(100000,1);
+ %! assert(min(x)>0); % *** Please report this!!! ***
+ %! assert(mean(x),1,0.01);
+ %! assert(var(x),1,0.03);
+ %! assert(skewness(x),2,0.06);
+ %! assert(kurtosis(x),6,0.7);
+ */
+
+ DEFUN_DLD (randg, args, ,
+ "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} randg (@var{a}, @var{x})\n\
+ @deftypefnx {Loadable Function} {} randg (@var{a}, @var{n}, @var{m})\n\
+ @deftypefnx {Loadable Function} {} randg (\"state\", @var{x})\n\
+ @deftypefnx {Loadable Function} {} randg (\"seed\", @var{x})\n\
+ Return a matrix with @code{gamma(@var{a},1)} distributed random elements.\n\
+ The arguments are handled the same as the arguments for @code{rand},\n\
+ except for the argument @var{a}.\n\
+ \n\
+ This can be used to generate many distributions:\n\
+ \n\
+ @table @asis\n\
+ @item @code{gamma (a,b)} for @code{a > -1}, @code{b > 0}\n\
+ @example\n\
+ r = b*randg(a)\n\
+ @end example\n\
+ @item @code{beta(a,b)} for @code{a > -1}, @code{b > -1}\n\
+ @example\n\
+ r1 = randg(a,1)\n\
+ r = r1 / (r1 + randg(b,1))\n\
+ @end example\n\
+ @item @code{Erlang(a, n)}\n\
+ @example\n\
+ r = a*randg(n)\n\
+ @end example\n\
+ @item @code{chisq(df)} for @code{df > 0}\n\
+ @example\n\
+ r = 2*randg(df/2)\n\
+ @end example\n\
+ @item @code{t(df)} for @code{0 < df < inf} (use randn if df is infinite)\n\
+ @example\n\
+ r = randn() / sqrt(2*randg(df/2)/df)\n\
+ @end example\n\
+ @item @code{F(n1,n2)} for @code{0 < n1}, @code{0 < n2}\n\
+ @example\n\
+ r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite\n\
+ r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite\n\
+ r = r1 / r2\n\n\
+ @end example\n\
+ @item negative @code{binomial (n, p)} for @code{n > 0}, @code{0 < p <= 1}\n\
+ @example\n\
+ r = randp((1-p)/p * randg(n))\n\
+ @end example\n\
+ @item non-central @code{chisq(df,L)}, for @code{df >= 0} and @code{L > 0}\n\
+ (use chisq if @code{L = 0})\n\
+ @example\n\
+ r = randp(L/2)\n\
+ r(r > 0) = 2*randg(r(r > 0))\n\
+ r(df > 0) += 2*randg(df(df > 0)/2)\n\
+ @end example\n\
+ @item @code{Dirichlet(a1,...,ak)}\n\
+ @example\n\
+ r = (randg(a1),...,randg(ak))\n\
+ r = r / sum(r)\n\
+ @end example\n\
+ @end table\n\
+ @seealso{rand,randn,rande,randp}\n\
+ @end deftypefn")
+ {
+ octave_value retval;
+
+ int nargin = args.length ();
+
+ if (nargin < 1)
+ error ("randg: insufficient arguments");
+ else
+ {
+ unwind_protect::begin_frame ("randg");
+
+ // This relies on the fact that elements are popped from the unwind
+ // stack in the reverse of the order they are pushed
+ // (i.e. current_distribution will be reset before calling
+ // reset_rand_generator()).
+
+ unwind_protect::add (reset_rand_generator, 0);
+ unwind_protect_str (current_distribution);
+
+ current_distribution = "gamma";
+
+ octave_rand::distribution (current_distribution);
+
+ retval = do_rand (args, nargin, "randg", true);
+
+ unwind_protect::run_frame ("randg");
+ }
+
+ return retval;
+ }
+
+ /*
+ %!test
+ %! rand("state",12)
+ %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! a=0.1; x = randg(a,100000,1);
+ %! assert(mean(x), a, 0.01);
+ %! assert(var(x), a, 0.01);
+ %! assert(skewness(x),2/sqrt(a), 1.);
+ %! assert(kurtosis(x),6/a, 50.);
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! a=0.95; x = randg(a,100000,1);
+ %! assert(mean(x), a, 0.01);
+ %! assert(var(x), a, 0.04);
+ %! assert(skewness(x),2/sqrt(a), 0.2);
+ %! assert(kurtosis(x),6/a, 2.);
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! a=1; x = randg(a,100000,1);
+ %! assert(mean(x),a, 0.01);
+ %! assert(var(x),a, 0.04);
+ %! assert(skewness(x),2/sqrt(a), 0.2);
+ %! assert(kurtosis(x),6/a, 2.);
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! a=10; x = randg(a,100000,1);
+ %! assert(mean(x), a, 0.1);
+ %! assert(var(x), a, 0.5);
+ %! assert(skewness(x),2/sqrt(a), 0.1);
+ %! assert(kurtosis(x),6/a, 0.5);
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! a=100; x = randg(a,100000,1);
+ %! assert(mean(x), a, 0.2);
+ %! assert(var(x), a, 2.);
+ %! assert(skewness(x),2/sqrt(a), 0.05);
+ %! assert(kurtosis(x),6/a, 0.2);
+ %!test
+ %! rand("seed",12)
+ %!assert(randg([-inf,-1,0,inf,nan]),[nan,nan,nan,nan,nan]) % *** Please report
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! a=0.1; x = randg(a,100000,1);
+ %! assert(mean(x), a, 0.01);
+ %! assert(var(x), a, 0.01);
+ %! assert(skewness(x),2/sqrt(a), 1.);
+ %! assert(kurtosis(x),6/a, 50.);
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! a=0.95; x = randg(a,100000,1);
+ %! assert(mean(x), a, 0.01);
+ %! assert(var(x), a, 0.04);
+ %! assert(skewness(x),2/sqrt(a), 0.2);
+ %! assert(kurtosis(x),6/a, 2.);
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! a=1; x = randg(a,100000,1);
+ %! assert(mean(x),a, 0.01);
+ %! assert(var(x),a, 0.04);
+ %! assert(skewness(x),2/sqrt(a), 0.2);
+ %! assert(kurtosis(x),6/a, 2.);
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! a=10; x = randg(a,100000,1);
+ %! assert(mean(x), a, 0.1);
+ %! assert(var(x), a, 0.5);
+ %! assert(skewness(x),2/sqrt(a), 0.1);
+ %! assert(kurtosis(x),6/a, 0.5);
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! a=100; x = randg(a,100000,1);
+ %! assert(mean(x), a, 0.2);
+ %! assert(var(x), a, 2.);
+ %! assert(skewness(x),2/sqrt(a), 0.05);
+ %! assert(kurtosis(x),6/a, 0.2);
+ */
+
+
+ DEFUN_DLD (randp, args, ,
+ "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} randp (@var{l}, @var{x})\n\
+ @deftypefnx {Loadable Function} {} randp (@var{l}, @var{n}, @var{m})\n\
+ @deftypefnx {Loadable Function} {} randp (\"state\", @var{x})\n\
+ @deftypefnx {Loadable Function} {} randp (\"seed\", @var{x})\n\
+ Return a matrix with Poisson distributed random elements. The arguments\n\
+ are handled the same as the arguments for @code{rand}, except for the\n\
+ argument @var{l}.\n\
+ \n\
+ Five different algorithms are used depending on the range of @var{l}\n\
+ and whether or not @var{l} is a scalar or a matrix.\n\
+ \n\
+ @table @asis\n\
+ @item For scalar @var{l} <= 12, use direct method.\n\
+ Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\
+ @item For scalar @var{l} > 12, use rejection method.[1]\n\
+ Press, et al., 'Numerical Recipes in C', Cambridge University Press, 1992.\n\
+ @item For matrix @var{l} <= 10, use inversion method.[2]\n\
+ Stadlober E., et al., WinRand source code, available via FTP.\n\
+ @item For matrix @var{l} > 10, use patchwork rejection method.\n\
+ Stadlober E., et al., WinRand source code, available via FTP, or\n\
+ H. Zechner, 'Efficient sampling from continuous and discrete\n\
+ unimodal distributions', Doctoral Dissertaion, 156pp., Technical\n\
+ University Graz, Austria, 1994.\n\
+ @item For @var{l} > 1e8, use normal approximation.\n\
+ L. Montanet, et al., 'Review of Particle Properties', Physical Review\n\
+ D 50 p1284, 1994\n\
+ @end table\n\
+ @seealso{rand,randn,rande,randg}\n\
+ @end deftypefn")
+ {
+ octave_value retval;
+
+ int nargin = args.length ();
+
+ if (nargin < 1)
+ error ("randp: insufficient arguments");
+ else
+ {
+ unwind_protect::begin_frame ("randp");
+
+ // This relies on the fact that elements are popped from the unwind
+ // stack in the reverse of the order they are pushed
+ // (i.e. current_distribution will be reset before calling
+ // reset_rand_generator()).
+
+ unwind_protect::add (reset_rand_generator, 0);
+ unwind_protect_str (current_distribution);
+
+ current_distribution = "poisson";
+
+ octave_rand::distribution (current_distribution);
+
+ retval = do_rand (args, nargin, "randp", true);
+
+ unwind_protect::run_frame ("randp");
+ }
+
+ return retval;
+ }
+
+ /*
+ %!test
+ %! rand("state",12)
+ %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! for a=[5 15]
+ %! x = randp(a,100000,1);
+ %! assert(min(x)>=0); % *** Please report this!!! ***
+ %! assert(mean(x),a,0.03);
+ %! assert(var(x),a,0.2);
+ %! assert(skewness(x),1/sqrt(a),0.03);
+ %! assert(kurtosis(x),1/a,0.08);
+ %! end
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! for a=[5 15]
+ %! x = randp(a*ones(100000,1),100000,1);
+ %! assert(min(x)>=0); % *** Please report this!!! ***
+ %! assert(mean(x),a,0.03);
+ %! assert(var(x),a,0.2);
+ %! assert(skewness(x),1/sqrt(a),0.03);
+ %! assert(kurtosis(x),1/a,0.08);
+ %! end
+ %!test
+ %! rand("seed",12)
+ %!assert(randp([-inf,-1,0,inf,nan]),[nan,nan,0,nan,nan]); % *** Please report
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! for a=[5 15]
+ %! x = randp(a,100000,1);
+ %! assert(min(x)>=0); % *** Please report this!!! ***
+ %! assert(mean(x),a,0.03);
+ %! assert(var(x),a,0.2);
+ %! assert(skewness(x),1/sqrt(a),0.03);
+ %! assert(kurtosis(x),1/a,0.08);
+ %! end
+ %!test
+ %! % statistical tests may fail occasionally.
+ %! for a=[5 15]
+ %! x = randp(a*ones(100000,1),100000,1);
+ %! assert(min(x)>=0); % *** Please report this!!! ***
+ %! assert(mean(x),a,0.03);
+ %! assert(var(x),a,0.2);
+ %! assert(skewness(x),1/sqrt(a),0.03);
+ %! assert(kurtosis(x),1/a,0.08);
+ %! end
+ */
+
+ /*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*** ./liboctave/randpoisson.c.orig5 2006-04-01 16:16:20.327039256 +0200
--- ./liboctave/randpoisson.c 2006-03-30 22:35:52.096397840 +0200
***************
*** 0 ****
--- 1,446 ----
+ /* This code is in the public domain */
+
+ /* Needs the following defines:
+ * NAN: value to return for Not-A-Number
+ * RUNI: uniform generator on (0,1)
+ * RNOR: normal generator
+ * LGAMMA: log gamma function
+ * INFINITE: function to test whether a value is infinite
+ */
+
+ #if defined (HAVE_CONFIG_H)
+ #include <config.h>
+ #endif
+
+ #include <math.h>
+ #include <stdio.h>
+
+ #include "f77-fcn.h"
+ #include "lo-ieee.h"
+ #include "lo-error.h"
+ #include "randmtzig.h"
+ #include "randpoisson.h"
+
+ #undef NAN
+ #define NAN octave_NaN
+ #define INFINITE lo_ieee_isinf
+ #define RUNI oct_randu()
+ #define RNOR oct_randn()
+ #define LGAMMA xlgamma
+
+ F77_RET_T
+ F77_FUNC (dlgams, DLGAMS) (const double *, double *, double *);
+
+ static double
+ xlgamma (double x)
+ {
+ double result;
+ double sgngam;
+
+ if (lo_ieee_isnan (x))
+ result = x;
+ else if (x <= 0 || lo_ieee_isinf (x))
+ result = octave_Inf;
+ else
+ F77_XFCN (dlgams, DLGAMS, (&x, &result, &sgngam));
+
+ return result;
+ }
+
+ /* ---- pprsc.c from Stadloeber's winrand --- */
+
+ #include <math.h>
+
+ /* flogfak(k) = ln(k!) */
+ static double
+ flogfak (double k)
+ {
+ #define C0 9.18938533204672742e-01
+ #define C1 8.33333333333333333e-02
+ #define C3 -2.77777777777777778e-03
+ #define C5 7.93650793650793651e-04
+ #define C7 -5.95238095238095238e-04
+
+ static double logfak[30L] = {
+ 0.00000000000000000, 0.00000000000000000, 0.69314718055994531,
+ 1.79175946922805500, 3.17805383034794562, 4.78749174278204599,
+ 6.57925121201010100, 8.52516136106541430, 10.60460290274525023,
+ 12.80182748008146961, 15.10441257307551530, 17.50230784587388584,
+ 19.98721449566188615, 22.55216385312342289, 25.19122118273868150,
+ 27.89927138384089157, 30.67186010608067280, 33.50507345013688888,
+ 36.39544520803305358, 39.33988418719949404, 42.33561646075348503,
+ 45.38013889847690803, 48.47118135183522388, 51.60667556776437357,
+ 54.78472939811231919, 58.00360522298051994, 61.26170176100200198,
+ 64.55753862700633106, 67.88974313718153498, 71.25703896716800901
+ };
+
+ double r, rr;
+
+ if (k >= 30.0)
+ {
+ r = 1.0 / k;
+ rr = r * r;
+ return ((k + 0.5)*log(k) - k + C0 + r*(C1 + rr*(C3 + rr*(C5 + rr*C7))));
+ }
+ else
+ return (logfak[(int)k]);
+ }
+
+
+ /******************************************************************
+ * *
+ * Poisson Distribution - Patchwork Rejection/Inversion *
+ * *
+ ******************************************************************
+ * *
+ * For parameter my < 10 Tabulated Inversion is applied. *
+ * For my >= 10 Patchwork Rejection is employed: *
+ * The area below the histogram function f(x) is rearranged in *
+ * its body by certain point reflections. Within a large center *
+ * interval variates are sampled efficiently by rejection from *
+ * uniform hats. Rectangular immediate acceptance regions speed *
+ * up the generation. The remaining tails are covered by *
+ * exponential functions. *
+ * *
+ ******************************************************************
+ * *
+ * FUNCTION : - pprsc samples a random number from the Poisson *
+ * distribution with parameter my > 0. *
+ * REFERENCE : - H. Zechner (1994): Efficient sampling from *
+ * continuous and discrete unimodal distributions, *
+ * Doctoral Dissertation, 156 pp., Technical *
+ * University Graz, Austria. *
+ * SUBPROGRAM : - drand(seed) ... (0,1)-Uniform generator with *
+ * unsigned long integer *seed. *
+ * *
+ * Implemented by H. Zechner, January 1994 *
+ * Revised by F. Niederl, July 1994 *
+ * *
+ ******************************************************************/
+
+ static double
+ f (double k, double l_nu, double c_pm)
+ {
+ return exp(k * l_nu - flogfak(k) - c_pm);
+ }
+
+ static double
+ pprsc (double my)
+ {
+ static double my_last = -1.0;
+ static double m, k2, k4, k1, k5;
+ static double dl, dr, r1, r2, r4, r5, ll, lr, l_my, c_pm,
+ f1, f2, f4, f5, p1, p2, p3, p4, p5, p6;
+ double Dk, X, Y;
+ double Ds, U, V, W;
+
+ if (my != my_last)
+ { /* set-up */
+ my_last = my;
+ /* approximate deviation of reflection points k2, k4 from my - 1/2 */
+ Ds = sqrt(my + 0.25);
+
+ /* mode m, reflection points k2 and k4, and points k1 and k5, */
+ /* which delimit the centre region of h(x) */
+ m = floor(my);
+ k2 = ceil(my - 0.5 - Ds);
+ k4 = floor(my - 0.5 + Ds);
+ k1 = k2 + k2 - m + 1L;
+ k5 = k4 + k4 - m;
+
+ /* range width of the critical left and right centre region */
+ dl = (k2 - k1);
+ dr = (k5 - k4);
+
+ /* recurrence constants r(k)=p(k)/p(k-1) at k = k1, k2, k4+1, k5+1 */
+ r1 = my / k1;
+ r2 = my / k2;
+ r4 = my / (k4 + 1.0);
+ r5 = my / (k5 + 1.0);
+
+ /* reciprocal values of the scale parameters of exp. tail envelope */
+ ll = log(r1); /* expon. tail left */
+ lr = -log(r5); /* expon. tail right*/
+
+ /* Poisson constants, necessary for computing function values f(k) */
+ l_my = log(my);
+ c_pm = m * l_my - flogfak(m);
+
+ /* function values f(k) = p(k)/p(m) at k = k2, k4, k1, k5 */
+ f2 = f(k2, l_my, c_pm);
+ f4 = f(k4, l_my, c_pm);
+ f1 = f(k1, l_my, c_pm);
+ f5 = f(k5, l_my, c_pm);
+
+ /* area of the two centre and the two exponential tail regions */
+ /* area of the two immediate acceptance regions between k2, k4 */
+ p1 = f2 * (dl + 1.0); /* immed. left */
+ p2 = f2 * dl + p1; /* centre left */
+ p3 = f4 * (dr + 1.0) + p2; /* immed. right */
+ p4 = f4 * dr + p3; /* centre right */
+ p5 = f1 / ll + p4; /* exp. tail left */
+ p6 = f5 / lr + p5; /* exp. tail right*/
+ }
+
+ for (;;)
+ {
+ /* generate uniform number U -- U(0, p6) */
+ /* case distinction corresponding to U */
+ if ((U = RUNI * p6) < p2)
+ { /* centre left */
+
+ /* immediate acceptance region
+ R2 = [k2, m) *[0, f2), X = k2, ... m -1 */
+ if ((V = U - p1) < 0.0) return(k2 + floor(U/f2));
+ /* immediate acceptance region
+ R1 = [k1, k2)*[0, f1), X = k1, ... k2-1 */
+ if ((W = V / dl) < f1 ) return(k1 + floor(V/f1));
+
+ /* computation of candidate X < k2, and its counterpart Y > k2 */
+ /* either squeeze-acceptance of X or acceptance-rejection of Y */
+ Dk = floor(dl * RUNI) + 1.0;
+ if (W <= f2 - Dk * (f2 - f2/r2))
+ { /* quick accept of */
+ return(k2 - Dk); /* X = k2 - Dk */
+ }
+ if ((V = f2 + f2 - W) < 1.0)
+ { /* quick reject of Y*/
+ Y = k2 + Dk;
+ if (V <= f2 + Dk * (1.0 - f2)/(dl + 1.0))
+ { /* quick accept of */
+ return(Y); /* Y = k2 + Dk */
+ }
+ if (V <= f(Y, l_my, c_pm)) return(Y); /* final accept of Y*/
+ }
+ X = k2 - Dk;
+ }
+ else if (U < p4)
+ { /* centre right */
+ /* immediate acceptance region
+ R3 = [m, k4+1)*[0, f4), X = m, ... k4 */
+ if ((V = U - p3) < 0.0) return(k4 - floor((U - p2)/f4));
+ /* immediate acceptance region
+ R4 = [k4+1, k5+1)*[0, f5) */
+ if ((W = V / dr) < f5 ) return(k5 - floor(V/f5));
+
+ /* computation of candidate X > k4, and its counterpart Y < k4 */
+ /* either squeeze-acceptance of X or acceptance-rejection of Y */
+ Dk = floor(dr * RUNI) + 1.0;
+ if (W <= f4 - Dk * (f4 - f4*r4))
+ { /* quick accept of */
+ return(k4 + Dk); /* X = k4 + Dk */
+ }
+ if ((V = f4 + f4 - W) < 1.0)
+ { /* quick reject of Y*/
+ Y = k4 - Dk;
+ if (V <= f4 + Dk * (1.0 - f4)/ dr)
+ { /* quick accept of */
+ return(Y); /* Y = k4 - Dk */
+ }
+ if (V <= f(Y, l_my, c_pm)) return(Y); /* final accept of Y*/
+ }
+ X = k4 + Dk;
+ }
+ else
+ {
+ W = RUNI;
+ if (U < p5)
+ { /* expon. tail left */
+ Dk = floor(1.0 - log(W)/ll);
+ if ((X = k1 - Dk) < 0L) continue; /* 0 <= X <= k1 - 1 */
+ W *= (U - p4) * ll; /* W -- U(0, h(x)) */
+ if (W <= f1 - Dk * (f1 - f1/r1))
+ return(X); /* quick accept of X*/
+ }
+ else
+ { /* expon. tail right*/
+ Dk = floor(1.0 - log(W)/lr);
+ X = k5 + Dk; /* X >= k5 + 1 */
+ W *= (U - p5) * lr; /* W -- U(0, h(x)) */
+ if (W <= f5 - Dk * (f5 - f5*r5))
+ return(X); /* quick accept of X*/
+ }
+ }
+
+ /* acceptance-rejection test of candidate X from the original area */
+ /* test, whether W <= f(k), with W = U*h(x) and U -- U(0, 1)*/
+ /* log f(X) = (X - m)*log(my) - log X! + log m! */
+ if (log(W) <= X * l_my - flogfak(X) - c_pm) return(X);
+ }
+ }
+ /* ---- pprsc.c end ------ */
+
+
+ /* The remainder of the file is by Paul Kienzle */
+
+ /* Given uniform u, find x such that CDF(L,x)==u. Return x. */
+ static void
+ poisson_cdf_lookup(double lambda, double *p, size_t n)
+ {
+ /* Table size is predicated on the maximum value of lambda
+ * we want to store in the table, and the maximum value of
+ * returned by the uniform random number generator on [0,1).
+ * With lambda==10 and u_max = 1 - 1/(2^32+1), we
+ * have poisson_pdf(lambda,36) < 1-u_max. If instead our
+ * generator uses more bits of mantissa or returns a value
+ * in the range [0,1], then for lambda==10 we need a table
+ * size of 46 instead. For long doubles, the table size
+ * will need to be longer still. */
+ #define TABLESIZE 46
+ double t[TABLESIZE];
+
+ /* Precompute the table for the u up to and including 0.458.
+ * We will almost certainly need it. */
+ int intlambda = (int)floor(lambda);
+ double P;
+ int tableidx;
+ size_t i = n;
+
+ t[0] = P = exp(-lambda);
+ for (tableidx = 1; tableidx <= intlambda; tableidx++) {
+ P = P*lambda/(double)tableidx;
+ t[tableidx] = t[tableidx-1] + P;
+ }
+
+ while (i-- > 0) {
+ double u = RUNI;
+
+ /* If u > 0.458 we know we can jump to floor(lambda) before
+ * comparing (this observation is based on Stadlober's winrand
+ * code). For lambda >= 1, this will be a win. Lambda < 1
+ * is already fast, so adding an extra comparison is not a
+ * problem. */
+ int k = (u > 0.458 ? intlambda : 0);
+
+ /* We aren't using a for loop here because when we find the
+ * right k we want to jump to the next iteration of the
+ * outer loop, and the continue statement will only work for
+ * the inner loop. */
+ nextk:
+ if ( u <= t[k] ) {
+ p[i] = (double) k;
+ continue;
+ }
+ if (++k < tableidx) goto nextk;
+
+ /* We only need high values of the table very rarely so we
+ * don't automatically compute the entire table. */
+ while (tableidx < TABLESIZE) {
+ P = P*lambda/(double)tableidx;
+ t[tableidx] = t[tableidx-1] + P;
+ /* Make sure we converge to 1.0 just in case u is uniform
+ * on [0,1] rather than [0,1). */
+ if (t[tableidx] == t[tableidx-1]) t[tableidx] = 1.0;
+ tableidx++;
+ if (u <= t[tableidx-1]) break;
+ }
+
+ /* We are assuming that the table size is big enough here.
+ * This should be true even if RUNI is returning values in
+ * the range [0,1] rather than [0,1).
+ */
+ p[i] = (double)(tableidx-1);
+ }
+ }
+
+ /* From Press, et al., Numerical Recipes */
+ static void
+ poisson_rejection (double lambda, double *p, size_t n)
+ {
+ double sq = sqrt(2.0*lambda);
+ double alxm = log(lambda);
+ double g = lambda*alxm - LGAMMA(lambda+1.0);
+ size_t i;
+
+ for (i = 0; i < n; i++)
+ {
+ double y, em, t;
+ do {
+ do {
+ y = tan(M_PI*RUNI);
+ em = sq * y + lambda;
+ } while (em < 0.0);
+ em = floor(em);
+ t = 0.9*(1.0+y*y)*exp(em*alxm-flogfak(em)-g);
+ } while (RUNI > t);
+ p[i] = em;
+ }
+ }
+
+ /* The cutoff of L <= 1e8 in the following two functions before using
+ * the normal approximation is based on:
+ * > L=1e8; x=floor(linspace(0,2*L,1000));
+ * > max(abs(normal_pdf(x,L,L)-poisson_pdf(x,L)))
+ * ans = 1.1376e-28
+ * For L=1e7, the max is around 1e-9, which is within the step size of RUNI.
+ * For L>1e10 the pprsc function breaks down, as I saw from the histogram
+ * of a large sample, so 1e8 is both small enough and large enough. */
+
+ /* Generate a set of poisson numbers with the same distribution */
+ void
+ oct_fill_randp (double L, octave_idx_type n, double *p)
+ {
+ octave_idx_type i;
+ if (L < 0.0 || INFINITE(L))
+ {
+ for (i=0; i<n; i++)
+ p[i] = NAN;
+ }
+ else if (L <= 10.0)
+ {
+ poisson_cdf_lookup(L, p, n);
+ }
+ else if (L <= 1e8)
+ {
+ for (i=0; i<n; i++)
+ p[i] = pprsc(L);
+ }
+ else
+ {
+ /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
+ const double sqrtL = sqrt(L);
+ for (i = 0; i < L; i++)
+ {
+ p[i] = floor(RNOR*sqrtL + L + 0.5);
+ if (p[i] < 0.0)
+ p[i] = 0.0; /* will probably never happen */
+ }
+ }
+ }
+
+ /* Generate one poisson variate */
+ double
+ oct_randp (double L)
+ {
+ double ret;
+ if (L < 0.0) ret = NAN;
+ else if (L <= 12.0) {
+ /* From Press, et al. Numerical recipes */
+ double g = exp(-L);
+ int em = -1;
+ double t = 1.0;
+ do {
+ ++em;
+ t *= RUNI;
+ } while (t > g);
+ ret = em;
+ } else if (L <= 1e8) {
+ /* numerical recipes */
+ poisson_rejection(L, &ret, 1);
+ } else if (INFINITE(L)) {
+ /* XXX FIXME XXX R uses NaN, but the normal approx. suggests that as
+ * limit should be inf. Which is correct? */
+ ret = NAN;
+ } else {
+ /* normal approximation: from Phys. Rev. D (1994) v50 p1284 */
+ ret = floor(RNOR*sqrt(L) + L + 0.5);
+ if (ret < 0.0) ret = 0.0; /* will probably never happen */
+ }
+ return ret;
+ }
+
+ /*
+ ;;; Local Variables: ***
+ ;;; mode: C ***
+ ;;; End: ***
+ */
*** ./liboctave/randpoisson.h.orig5 2006-04-01 16:16:22.443717472 +0200
--- ./liboctave/randpoisson.h 2006-03-30 22:34:43.295857112 +0200
***************
*** 0 ****
--- 1,24 ----
+ /* This code is in the public domain */
+
+ #ifndef _RANDPOISSON_H
+
+ #include "oct-types.h"
+
+ #ifdef __cplusplus
+ extern "C" {
+ #endif
+
+ extern double oct_randp (double L);
+ extern void oct_fill_randp (double L, octave_idx_type n, double *p);
+
+ #ifdef __cplusplus
+ }
+ #endif
+ #endif
+
+ /*
+ ;;; Local Variables: ***
+ ;;; mode: C ***
+ ;;; End: ***
+ */
+
*** ./liboctave/randgamma.c.orig5 2006-04-01 16:16:31.251378504 +0200
--- ./liboctave/randgamma.c 2006-03-30 22:37:16.879508848 +0200
***************
*** 0 ****
--- 1,126 ----
+ /* This code is in the public domain */
+
+ /*
+
+ double randg(a)
+ void fill_randg(a,n,x)
+
+ Generate a series of standard gamma distributions.
+
+ See: Marsaglia G and Tsang W (2000), "A simple method for generating
+ gamma variables", ACM Transactions on Mathematical Software 26(3) 363-372
+
+ Needs the following defines:
+ * NAN: value to return for Not-A-Number
+ * RUNI: uniform generator on (0,1)
+ * RNOR: normal generator
+ * REXP: exponential generator, or -log(RUNI) if one isn't available
+ * INFINITE: function to test whether a value is infinite
+
+ Test using:
+ mean = a
+ variance = a
+ skewness = 2/sqrt(a)
+ kurtosis = 3 + 6/sqrt(a)
+
+ Note that randg can be used to generate many distributions:
+
+ gamma(a,b) for a>0, b>0 (from R)
+ r = b*randg(a)
+ beta(a,b) for a>0, b>0
+ r1 = randg(a,1)
+ r = r1 / (r1 + randg(b,1))
+ Erlang(a,n)
+ r = a*randg(n)
+ chisq(df) for df>0
+ r = 2*randg(df/2)
+ t(df) for 0<df<inf (use randn if df is infinite)
+ r = randn() / sqrt(2*randg(df/2)/df)
+ F(n1,n2) for 0<n1, 0<n2
+ r1 = 2*randg(n1/2)/n1 or 1 if n1 is infinite
+ r2 = 2*randg(n2/2)/n2 or 1 if n2 is infinite
+ r = r1 / r2
+ negative binonial (n, p) for n>0, 0<p<=1
+ r = randp((1-p)/p * randg(n))
+ (from R, citing Devroye(1986), Non-Uniform Random Variate Generation)
+ non-central chisq(df,L), for df>=0 and L>0 (use chisq if L=0)
+ r = randp(L/2)
+ r(r>0) = 2*randg(r(r>0))
+ r(df>0) += 2*randg(df(df>0)/2)
+ (from R, citing formula 29.5b-c in Johnson, Kotz, Balkrishnan(1995))
+ Dirichlet(a1,...,ak) for ai > 0
+ r = (randg(a1),...,randg(ak))
+ r = r / sum(r)
+ (from GSL, citing Law & Kelton(1991), Simulation Modeling and Analysis)
+ */
+
+ #if defined (HAVE_CONFIG_H)
+ #include <config.h>
+ #endif
+
+ #include <math.h>
+ #include <stdio.h>
+
+ #include "lo-ieee.h"
+ #include "randmtzig.h"
+ #include "randgamma.h"
+
+ #undef NAN
+ #define NAN octave_NaN
+ #define INFINITE lo_ieee_isinf
+ #define RUNI oct_randu()
+ #define RNOR oct_randn()
+ #define REXP oct_rande()
+
+ void
+ oct_fill_randg (double a, octave_idx_type n, double *r)
+ {
+ octave_idx_type i;
+ /* If a < 1, start by generating gamma(1+a) */
+ const double d = (a < 1. ? 1.+a : a) - 1./3.;
+ const double c = 1./sqrt(9.*d);
+
+ /* Handle invalid cases */
+ if (a <= 0 || INFINITE(a))
+ {
+ for (i=0; i < n; i++)
+ r[i] = NAN;
+ return;
+ }
+
+ for (i=0; i < n; i++)
+ {
+ double x, xsq, v, u;
+ restart:
+ x = RNOR;
+ v = (1+c*x);
+ v *= v*v;
+ if (v <= 0)
+ goto restart; /* rare, so don't bother moving up */
+ u = RUNI;
+ xsq = x*x;
+ if (u >= 1.-0.0331*xsq*xsq && log(u) >= 0.5*xsq + d*(1-v+log(v)))
+ goto restart;
+ r[i] = d*v;
+ }
+ if (a < 1)
+ { /* Use gamma(a) = gamma(1+a)*U^(1/a) */
+ /* Given REXP = -log(U) then U^(1/a) = exp(-REXP/a) */
+ for (i = 0; i < n; i++)
+ r[i] *= exp(-REXP/a);
+ }
+ }
+
+ double
+ oct_randg (double a)
+ {
+ double ret;
+ oct_fill_randg(a,1,&ret);
+ return ret;
+ }
+
+ /*
+ ;;; Local Variables: ***
+ ;;; mode: C ***
+ ;;; End: ***
+ */
*** ./liboctave/randgamma.h.orig5 2006-04-01 16:16:40.124029656 +0200
--- ./liboctave/randgamma.h 2006-03-30 22:36:25.275353872 +0200
***************
*** 0 ****
--- 1,24 ----
+ /* This code is in the public domain */
+
+ #ifndef _RANDGAMMA_H
+
+ #include "oct-types.h"
+
+ #ifdef __cplusplus
+ extern "C" {
+ #endif
+
+ extern double oct_randg (double a);
+ extern void oct_fill_randg (double a, octave_idx_type n, double *p);
+
+ #ifdef __cplusplus
+ }
+ #endif
+ #endif
+
+ /*
+ ;;; Local Variables: ***
+ ;;; mode: C ***
+ ;;; End: ***
+ */
+
*** ./liboctave/randmtzig.c.orig5 2006-04-01 16:16:46.463065976 +0200
--- ./liboctave/randmtzig.c 2006-03-30 23:51:16.630564096 +0200
***************
*** 0 ****
--- 1,686 ----
+ /*
+ A C-program for MT19937, with initialization improved 2002/2/10.
+ Coded by Takuji Nishimura and Makoto Matsumoto.
+ This is a faster version by taking Shawn Cokus's optimization,
+ Matthe Bellew's simplification, Isaku Wada's real version.
+ David Bateman added normal and exponential distributions following
+ Marsaglia and Tang's Ziggurat algorithm.
+
+ Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+ Copyright (C) 2004, David Bateman
+ All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+
+ 1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ 2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+
+ 3. The names of its contributors may not be used to endorse or promote
+ products derived from this software without specific prior written
+ permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER
+ OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+ Any feedback is very welcome.
+ http://www.math.keio.ac.jp/matumoto/emt.html
+ email: address@hidden
+
+ * 2006-04-01 David Bateman
+ * * convert for use in octave, declaring static functions only used
+ * here and adding oct_ to functions visible externally
+ * * inverse sense of ALLBITS
+ * 2004-01-19 Paul Kienzle
+ * * comment out main
+ * add init_by_entropy, get_state, set_state
+ * * converted to allow compiling by C++ compiler
+ *
+ * 2004-01-25 David Bateman
+ * * Add Marsaglia and Tsang Ziggurat code
+ *
+ * 2004-07-13 Paul Kienzle
+ * * make into an independent library with some docs.
+ * * introduce new main and test code.
+ *
+ * 2004-07-28 Paul Kienzle & David Bateman
+ * * add -DALLBITS flag for 32 vs. 53 bits of randomness in mantissa
+ * * make the naming scheme more uniform
+ * * add -DHAVE_X86 for faster support of 53 bit mantissa on x86 arch.
+ *
+ * 2005-02-23 Paul Kienzle
+ * * fix -DHAVE_X86_32 flag and add -DUSE_X86_32=0|1 for explicit control
+ */
+
+ /*
+ === Build instructions ===
+
+ Compile with -DHAVE_GETTIMEOFDAY if the gettimeofday function is
+ available. This is not necessary if your architecture has
+ /dev/urandom defined.
+
+ Compile with -DALLBITS to disable 53-bit random numbers. This is about
+ 50% slower than using 32-bit random numbers.
+
+ Uses implicit -Di386 or explicit -DHAVE_X86_32 to determine if CPU=x86.
+ You can force X86 behaviour with -DUSE_X86_32=1, or suppress it with
+ -DUSE_X86_32=0. You should also consider -march=i686 or similar for
+ extra performance. Check whether -DUSE_X86_32=0 is faster on 64-bit
+ x86 architectures.
+
+ If you want to replace the Mersenne Twister with another
+ generator then redefine randi32 appropriately.
+
+ === Usage instructions ===
+ Before using any of the generators, initialize the state with one of
+ oct_init_by_int, oct_init_by_array or oct_init_by_entropy.
+
+ All generators share the same state vector.
+
+ === Mersenne Twister ===
+ void oct_init_by_int(uint32_t s) 32-bit initial state
+ void oct_init_by_array(uint32_t k[],int m) m*32-bit initial state
+ void oct_init_by_entropy(void) random initial state
+ void oct_get_state(uint32_t save[MT_N+1]) saves state in array
+ void oct_set_state(uint32_t save[MT_N+1]) restores state from array
+ static inline uint32_t randmt(void) returns 32-bit unsigned int
+
+ === inline generators ===
+ static inline uint32_t randi32(void) returns 32-bit unsigned int
+ static inline uint64_t randi53(void) returns 53-bit unsigned int
+ static inline uint64_t randi54(void) returns 54-bit unsigned int
+ static inline uint64_t randi64(void) returns 64-bit unsigned int
+ static inline double randu32(void) returns 32-bit uniform in (0,1)
+ static inline double randu53(void) returns 53-bit uniform in (0,1)
+
+ double oct_randu(void) returns M-bit uniform in (0,1)
+ double oct_randn(void) returns M-bit standard normal
+ double oct_rande(void) returns N-bit standard exponential
+
+ === Array generators ===
+ void oct_fill_randi32(octave_idx_type, uint32_t [])
+ void oct_fill_randi64(octave_idx_type, uint64_t [])
+ void oct_fill_randu(octave_idx_type, double [])
+ void oct_fill_randn(octave_idx_type, double [])
+ void oct_fill_rande(octave_idx_type, double [])
+
+ */
+
+ #if defined (HAVE_CONFIG_H)
+ #include <config.h>
+ #endif
+
+ #include <math.h>
+ #include <stdio.h>
+ #include <time.h>
+
+ #ifdef HAVE_GETTIMEOFDAY
+ #include <sys/time.h>
+ #endif
+
+ #include "randmtzig.h"
+
+ /* XXX FIXME XXX may want to suppress X86 if sizeof(long)>4 */
+ #if !defined(USE_X86_32)
+ # if defined(i386) || defined(HAVE_X86_32)
+ # define USE_X86_32 1
+ # else
+ # define USE_X86_32 0
+ # endif
+ #endif
+
+ /* ===== Mersenne Twister 32-bit generator ===== */
+
+ #define MT_M 397
+ #define MATRIX_A 0x9908b0dfUL /* constant vector a */
+ #define UMASK 0x80000000UL /* most significant w-r bits */
+ #define LMASK 0x7fffffffUL /* least significant r bits */
+ #define MIXBITS(u,v) ( ((u) & UMASK) | ((v) & LMASK) )
+ #define TWIST(u,v) ((MIXBITS(u,v) >> 1) ^ ((v)&1UL ? MATRIX_A : 0UL))
+
+ static uint32_t *next;
+ static uint32_t state[MT_N]; /* the array for the state vector */
+ static int left = 1;
+ static int initf = 0;
+ static int initt = 1;
+
+ /* initializes state[MT_N] with a seed */
+ void
+ oct_init_by_int (uint32_t s)
+ {
+ int j;
+ state[0] = s & 0xffffffffUL;
+ for (j = 1; j < MT_N; j++) {
+ state[j] = (1812433253UL * (state[j-1] ^ (state[j-1] >> 30)) + j);
+ /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
+ /* In the previous versions, MSBs of the seed affect */
+ /* only MSBs of the array state[]. */
+ /* 2002/01/09 modified by Makoto Matsumoto */
+ state[j] &= 0xffffffffUL; /* for >32 bit machines */
+ }
+ left = 1;
+ initf = 1;
+ }
+
+ /* initialize by an array with array-length */
+ /* init_key is the array for initializing keys */
+ /* key_length is its length */
+ void
+ oct_init_by_array (uint32_t init_key[], int key_length)
+ {
+ int i, j, k;
+ oct_init_by_int (19650218UL);
+ i = 1;
+ j = 0;
+ k = (MT_N > key_length ? MT_N : key_length);
+ for (; k; k--)
+ {
+ state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) * 1664525UL))
+ + init_key[j] + j; /* non linear */
+ state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+ i++;
+ j++;
+ if (i >= MT_N)
+ {
+ state[0] = state[MT_N-1];
+ i = 1;
+ }
+ if (j >= key_length)
+ j = 0;
+ }
+ for (k = MT_N - 1; k; k--)
+ {
+ state[i] = (state[i] ^ ((state[i-1] ^ (state[i-1] >> 30)) *
1566083941UL))
+ - i; /* non linear */
+ state[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
+ i++;
+ if (i >= MT_N)
+ {
+ state[0] = state[MT_N-1];
+ i = 1;
+ }
+ }
+
+ state[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
+ left = 1;
+ initf = 1;
+ }
+
+ void
+ oct_init_by_entropy (void)
+ {
+ uint32_t entropy[MT_N];
+ int n = 0;
+
+ /* Look for entropy in /dev/urandom */
+ FILE* urandom =fopen("/dev/urandom", "rb");
+ if (urandom)
+ {
+ while (n < MT_N)
+ {
+ unsigned char word[4];
+ if (fread(word, 4, 1, urandom) != 1)
+ break;
+ entropy[n++] = word[0]+(word[1]<<8)+(word[2]<<16)+(word[3]<<24);
+ }
+ fclose(urandom);
+ }
+
+ /* If there isn't enough entropy, gather some from various sources */
+ if (n < MT_N)
+ entropy[n++] = time(NULL); /* Current time in seconds */
+ if (n < MT_N)
+ entropy[n++] = clock(); /* CPU time used (usec) */
+ #ifdef HAVE_GETTIMEOFDAY
+ if (n < MT_N)
+ {
+ struct timeval tv;
+ if (gettimeofday(&tv, NULL) != -1)
+ entropy[n++] = tv.tv_usec; /* Fractional part of current time */
+ }
+ #endif
+ /* Send all the entropy into the initial state vector */
+ oct_init_by_array(entropy,n);
+ }
+
+ void
+ oct_set_state (uint32_t save[])
+ {
+ int i;
+ for (i=0; i < MT_N; i++)
+ state[i] = save[i];
+ left = save[MT_N];
+ next = state + (MT_N - left + 1);
+ }
+
+ void
+ oct_get_state (uint32_t save[])
+ {
+ int i;
+ for (i = 0; i < MT_N; i++)
+ save[i] = state[i];
+ save[MT_N] = left;
+ }
+
+ static void
+ next_state (void)
+ {
+ uint32_t *p = state;
+ int j;
+
+ /* if init_by_int() has not been called, */
+ /* a default initial seed is used */
+ /* if (initf==0) init_by_int(5489UL); */
+ /* Or better yet, a random seed! */
+ if (initf == 0)
+ oct_init_by_entropy();
+
+ left = MT_N;
+ next = state;
+
+ for (j = MT_N - MT_M + 1; --j; p++)
+ *p = p[MT_M] ^ TWIST(p[0], p[1]);
+
+ for (j = MT_M; --j; p++)
+ *p = p[MT_M-MT_N] ^ TWIST(p[0], p[1]);
+
+ *p = p[MT_M-MT_N] ^ TWIST(p[0], state[0]);
+ }
+
+ /* generates a random number on [0,0xffffffff]-interval */
+ static inline uint32_t
+ randmt (void)
+ {
+ register uint32_t y;
+
+ if (--left == 0)
+ next_state();
+ y = *next++;
+
+ /* Tempering */
+ y ^= (y >> 11);
+ y ^= (y << 7) & 0x9d2c5680UL;
+ y ^= (y << 15) & 0xefc60000UL;
+ return (y ^ (y >> 18));
+ }
+
+ /* ===== Uniform generators ===== */
+
+ /* Select which 32 bit generator to use */
+ #define randi32 randmt
+
+ static inline uint64_t
+ randi53 (void)
+ {
+ const uint32_t lo = randi32();
+ const uint32_t hi = randi32()&0x1FFFFF;
+ #if HAVE_X86_32
+ uint64_t u;
+ uint32_t *p = (uint32_t *)&u;
+ p[0] = lo;
+ p[1] = hi;
+ return u;
+ #else
+ return (((uint64_t)hi<<32)|lo);
+ #endif
+ }
+
+ static inline uint64_t
+ randi54 (void)
+ {
+ const uint32_t lo = randi32();
+ const uint32_t hi = randi32()&0x3FFFFF;
+ #if HAVE_X86_32
+ uint64_t u;
+ uint32_t *p = (uint32_t *)&u;
+ p[0] = lo;
+ p[1] = hi;
+ return u;
+ #else
+ return (((uint64_t)hi<<32)|lo);
+ #endif
+ }
+
+ static inline uint64_t
+ randi64 (void)
+ {
+ const uint32_t lo = randi32();
+ const uint32_t hi = randi32();
+ #if HAVE_X86_32
+ uint64_t u;
+ uint32_t *p = (uint32_t *)&u;
+ p[0] = lo;
+ p[1] = hi;
+ return u;
+ #else
+ return (((uint64_t)hi<<32)|lo);
+ #endif
+ }
+
+ /* generates a random number on (0,1)-real-interval */
+ static inline double
+ randu32 (void)
+ {
+ return ((double)randi32() + 0.5) * (1.0/4294967296.0);
+ /* divided by 2^32 */
+ }
+
+ /* generates a random number on (0,1) with 53-bit resolution */
+ static inline double
+ randu53 (void)
+ {
+ const uint32_t a=randi32()>>5;
+ const uint32_t b=randi32()>>6;
+ return(a*67108864.0+b+0.4) * (1.0/9007199254740992.0);
+ }
+
+ /* Determine mantissa for uniform doubles */
+ #ifdef ALLBITS
+ double
+ oct_randu (void)
+ {
+ return randu32();
+ }
+ #else
+ double
+ oct_randu (void)
+ {
+ return randu53();
+ }
+ #endif
+
+ /* ===== Ziggurat normal and exponential generators ===== */
+ #ifdef ALLBITS
+ # define ZIGINT uint32_t
+ # define EMANTISSA 4294967296.0 /* 32 bit mantissa */
+ # define ERANDI randi32() /* 32 bits for mantissa */
+ # define NMANTISSA 2147483648.0 /* 31 bit mantissa */
+ # define NRANDI randi32() /* 31 bits for mantissa + 1 bit sign */
+ # define RANDU randu32()
+ #else
+ # define ZIGINT uint64_t
+ # define EMANTISSA 9007199254740992.0 /* 53 bit mantissa */
+ # define ERANDI randi53() /* 53 bits for mantissa */
+ # define NMANTISSA EMANTISSA
+ # define NRANDI randi54() /* 53 bits for mantissa + 1 bit sign */
+ # define RANDU randu53()
+ #endif
+
+ #define ZIGGURAT_TABLE_SIZE 256
+
+ #define ZIGGURAT_NOR_R 3.6541528853610088
+ #define ZIGGURAT_NOR_INV_R 0.27366123732975828
+ #define NOR_SECTION_AREA 0.00492867323399
+
+ #define ZIGGURAT_EXP_R 7.69711747013104972
+ #define ZIGGURAT_EXP_INV_R 0.129918765548341586
+ #define EXP_SECTION_AREA 0.0039496598225815571993
+
+ static ZIGINT ki[ZIGGURAT_TABLE_SIZE];
+ static double wi[ZIGGURAT_TABLE_SIZE], fi[ZIGGURAT_TABLE_SIZE];
+ static ZIGINT ke[ZIGGURAT_TABLE_SIZE];
+ static double we[ZIGGURAT_TABLE_SIZE], fe[ZIGGURAT_TABLE_SIZE];
+
+ /*
+ This code is based on the paper Marsaglia and Tsang, "The ziggurat method
+ for generating random variables", Journ. Statistical Software. Code was
+ presented in this paper for a Ziggurat of 127 levels and using a 32 bit
+ integer random number generator. This version of the code, uses the
+ Mersenne Twister as the integer generator and uses 256 levels in the
+ Ziggurat. This has several advantages.
+
+ 1) As Marsaglia and Tsang themselves states, the more levels the few
+ times the expensive tail algorithm must be called
+ 2) The cycle time of the generator is determined by the integer
+ generator, thus the use of a Mersenne Twister for the core random
+ generator makes this cycle extremely long.
+ 3) The license on the original code was unclear, thus rewriting the code
+ from the article means we are free of copyright issues.
+ 4) Compile flag for full 53-bit random mantissa.
+
+ It should be stated that the authors made my life easier, by the fact that
+ the algorithm developed in the text of the article is for a 256 level
+ ziggurat, even if the code itself isn't...
+
+ One modification to the algorithm developed in the article, is that it is
+ assumed that 0 <= x < Inf, and "unsigned long"s are used, thus resulting in
+ terms like 2^32 in the code. As the normal distribution is defined between
+ -Inf < x < Inf, we effectively only have 31 bit integers plus a sign. Thus
+ in Marsaglia and Tsang, terms like 2^32 become 2^31. We use NMANTISSA for
+ this term. The exponential distribution is one sided so we use the
+ full 32 bits. We use EMANTISSA for this term.
+
+ It appears that I'm slightly slower than the code in the article, this
+ is partially due to a better generator of random integers than they
+ use. But might also be that the case of rapid return was optimized by
+ inlining the relevant code with a #define. As the basic Mersenne
+ Twister is only 25% faster than this code I suspect that the main
+ reason is just the use of the Mersenne Twister and not the inlining,
+ so I'm not going to try and optimize further.
+ */
+
+ static void
+ create_ziggurat_tables (void)
+ {
+ int i;
+ double x, x1;
+
+ /* Ziggurat tables for the normal distribution */
+ x1 = ZIGGURAT_NOR_R;
+ wi[255] = x1 / NMANTISSA;
+ fi[255] = exp (-0.5 * x1 * x1);
+
+ /* Index zero is special for tail strip, where Marsaglia and Tsang
+ * defines this as
+ * k_0 = 2^31 * r * f(r) / v, w_0 = 0.5^31 * v / f(r), f_0 = 1,
+ * where v is the area of each strip of the ziggurat.
+ */
+ ki[0] = (ZIGINT) (x1 * fi[255] / NOR_SECTION_AREA * NMANTISSA);
+ wi[0] = NOR_SECTION_AREA / fi[255] / NMANTISSA;
+ fi[0] = 1.;
+
+ for (i = 254; i > 0; i--)
+ {
+ /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
+ * need inverse operator of y = exp(-0.5*x*x) -> x = sqrt(-2*ln(y))
+ */
+ x = sqrt(-2. * log(NOR_SECTION_AREA / x1 + fi[i+1]));
+ ki[i+1] = (ZIGINT)(x / x1 * NMANTISSA);
+ wi[i] = x / NMANTISSA;
+ fi[i] = exp (-0.5 * x * x);
+ x1 = x;
+ }
+
+ ki[1] = 0;
+
+ /* Zigurrat tables for the exponential distribution */
+ x1 = ZIGGURAT_EXP_R;
+ we[255] = x1 / EMANTISSA;
+ fe[255] = exp (-x1);
+
+ /* Index zero is special for tail strip, where Marsaglia and Tsang
+ * defines this as
+ * k_0 = 2^32 * r * f(r) / v, w_0 = 0.5^32 * v / f(r), f_0 = 1,
+ * where v is the area of each strip of the ziggurat.
+ */
+ ke[0] = (ZIGINT) (x1 * fe[255] / EXP_SECTION_AREA * EMANTISSA);
+ we[0] = EXP_SECTION_AREA / fe[255] / EMANTISSA;
+ fe[0] = 1.;
+
+ for (i = 254; i > 0; i--)
+ {
+ /* New x is given by x = f^{-1}(v/x_{i+1} + f(x_{i+1})), thus
+ * need inverse operator of y = exp(-x) -> x = -ln(y)
+ */
+ x = - log(EXP_SECTION_AREA / x1 + fe[i+1]);
+ ke[i+1] = (ZIGINT)(x / x1 * EMANTISSA);
+ we[i] = x / EMANTISSA;
+ fe[i] = exp (-x);
+ x1 = x;
+ }
+ ke[1] = 0;
+
+ initt = 0;
+ }
+
+ /*
+ * Here is the guts of the algorithm. As Marsaglia and Tsang state the
+ * algorithm in their paper
+ *
+ * 1) Calculate a random signed integer j and let i be the index
+ * provided by the rightmost 8-bits of j
+ * 2) Set x = j * w_i. If j < k_i return x
+ * 3) If i = 0, then return x from the tail
+ * 4) If [f(x_{i-1}) - f(x_i)] * U < f(x) - f(x_i), return x
+ * 5) goto step 1
+ *
+ * Where f is the functional form of the distribution, which for a normal
+ * distribution is exp(-0.5*x*x)
+ */
+
+ double
+ oct_randn (void)
+ {
+ if (initt)
+ create_ziggurat_tables();
+
+ while (1)
+ {
+ /* The following code is specialized for 32-bit mantissa.
+ * Compared to the arbitrary mantissa code, there is a performance
+ * gain for 32-bits: PPC: 2%, MIPS: 8%, x86: 40%
+ * There is a bigger performance gain compared to using a full
+ * 53-bit mantissa: PPC: 60%, MIPS: 65%, x86: 240%
+ * Of course, different compilers and operating systems may
+ * have something to do with this.
+ */
+ #if !defined(ALLBITS)
+ # if HAVE_X86_32
+ /* 53-bit mantissa, 1-bit sign, x86 32-bit architecture */
+ double x;
+ int si,idx;
+ register uint32_t lo, hi;
+ int64_t rabs;
+ uint32_t *p = (uint32_t *)&rabs;
+ lo = randi32();
+ idx = lo&0xFF;
+ hi = randi32();
+ si = hi&UMASK;
+ p[0] = lo;
+ p[1] = hi&0x1FFFFF;
+ x = ( si ? -rabs : rabs ) * wi[idx];
+ # else /* !HAVE_X86_32 */
+ /* arbitrary mantissa (selected by NRANDI, with 1 bit for sign) */
+ const uint64_t r = NRANDI;
+ const int64_t rabs=r>>1;
+ const int idx = (int)(rabs&0xFF);
+ const double x = ( r&1 ? -rabs : rabs) * wi[idx];
+ # endif /* !HAVE_X86_32 */
+ if (rabs < (int64_t)ki[idx])
+ #else /* ALLBITS */
+ /* 32-bit mantissa */
+ const uint32_t r = randi32();
+ const uint32_t rabs = r&LMASK;
+ const int idx = (int)(r&0xFF);
+ const double x = ((int32_t)r) * wi[idx];
+ if (rabs < ki[idx])
+ #endif /* ALLBITS */
+ return x; /* 99.3% of the time we return here 1st try */
+ else if (idx == 0)
+ {
+ /* As stated in Marsaglia and Tsang
+ *
+ * For the normal tail, the method of Marsaglia[5] provides:
+ * generate x = -ln(U_1)/r, y = -ln(U_2), until y+y > x*x,
+ * then return r+x. Except that r+x is always in the positive
+ * tail!!!! Any thing random might be used to determine the
+ * sign, but as we already have r we might as well use it
+ *
+ * [PAK] but not the bottom 8 bits, since they are all 0 here!
+ */
+ double xx, yy;
+ do
+ {
+ xx = - ZIGGURAT_NOR_INV_R * log (RANDU);
+ yy = - log (RANDU);
+ }
+ while ( yy+yy <= xx*xx);
+ return (rabs&0x100 ? -ZIGGURAT_NOR_R-xx : ZIGGURAT_NOR_R+xx);
+ }
+ else if ((fi[idx-1] - fi[idx]) * RANDU + fi[idx] < exp(-0.5*x*x))
+ return x;
+ }
+ }
+
+ double
+ oct_rande (void)
+ {
+ if (initt)
+ create_ziggurat_tables();
+
+ while (1)
+ {
+ ZIGINT ri = ERANDI;
+ const int idx = (int)(ri & 0xFF);
+ const double x = ri * we[idx];
+ if (ri < ke[idx])
+ return x; // 98.9% of the time we return here 1st try
+ else if (idx == 0)
+ {
+ /* As stated in Marsaglia and Tsang
+ *
+ * For the exponential tail, the method of Marsaglia[5] provides:
+ * x = r - ln(U);
+ */
+ return ZIGGURAT_EXP_R - log(RANDU);
+ }
+ else if ((fe[idx-1] - fe[idx]) * RANDU + fe[idx] < exp(-x))
+ return x;
+ }
+ }
+
+ /* Array generators */
+ void
+ oct_fill_randu (octave_idx_type n, double *p)
+ {
+ octave_idx_type i;
+ for (i = 0; i < n; i++)
+ p[i] = oct_randu();
+ }
+
+ void
+ oct_fill_randn (octave_idx_type n, double *p)
+ {
+ octave_idx_type i;
+ for (i = 0; i < n; i++)
+ p[i] = oct_randn();
+ }
+
+ void
+ oct_fill_rande (octave_idx_type n, double *p)
+ {
+ octave_idx_type i;
+ for (i = 0; i < n; i++)
+ p[i] = oct_rande();
+ }
+
+ /*
+ ;;; Local Variables: ***
+ ;;; mode: C ***
+ ;;; End: ***
+ */
*** ./liboctave/randmtzig.h.orig5 2006-04-01 16:16:49.160655880 +0200
--- ./liboctave/randmtzig.h 2006-03-30 22:32:21.183461480 +0200
***************
*** 0 ****
--- 1,104 ----
+ /*
+ A C-program for MT19937, with initialization improved 2002/2/10.
+ Coded by Takuji Nishimura and Makoto Matsumoto.
+ This is a faster version by taking Shawn Cokus's optimization,
+ Matthe Bellew's simplification, Isaku Wada's real version.
+ David Bateman added normal and exponential distributions following
+ Marsaglia and Tang's Ziggurat algorithm.
+
+ Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
+ Copyright (C) 2004, David Bateman
+ All rights reserved.
+
+ Redistribution and use in source and binary forms, with or without
+ modification, are permitted provided that the following conditions
+ are met:
+
+ 1. Redistributions of source code must retain the above copyright
+ notice, this list of conditions and the following disclaimer.
+
+ 2. Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+
+ 3. The names of its contributors may not be used to endorse or promote
+ products derived from this software without specific prior written
+ permission.
+
+ THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+ A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER
+ OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+ */
+
+ #ifndef _RANDMTZIG_H
+ #define _RANDMTZIG_H
+
+ #ifdef HAVE_CONFIG_H
+ #include <config.h>
+ #endif
+
+ #include "oct-types.h"
+
+ #define MT_N 624
+
+ #ifdef __cplusplus
+ extern "C" {
+ #endif
+
+ #ifdef HAVE_INTTYPES_H
+ #include <inttypes.h>
+ #else
+ #if SIZEOF_INT == 4
+ typedef unsigned int uint32_t;
+ #elif SIZEOF_LONG == 4
+ typedef unsigned long uint32_t;
+ #else
+ #error "No 4 byte integer type found!"
+ #endif
+
+ #if SIZEOF_LONG == 8
+ typedef unsigned long uint64_t;
+ #else
+ #if SIZEOF_LONG_LONG == 8
+ typedef unsigned long long uint64_t;
+ #endif
+ #endif
+ #endif
+
+ /* === Mersenne Twister === */
+ extern void oct_init_by_int (uint32_t s);
+ extern void oct_init_by_array (uint32_t init_key[], int key_length);
+ extern void oct_init_by_entropy (void);
+ extern void oct_set_state (uint32_t save[]);
+ extern void oct_get_state (uint32_t save[]);
+
+ /* === Array generators === */
+ extern double oct_randu (void);
+ extern double oct_randn (void);
+ extern double oct_rande (void);
+
+ /* === Array generators === */
+ extern void oct_fill_randu (octave_idx_type n, double *p);
+ extern void oct_fill_randn (octave_idx_type n, double *p);
+ extern void oct_fill_rande (octave_idx_type n, double *p);
+
+ #ifdef __cplusplus
+ }
+ #endif
+ #endif
+
+ /*
+ ;;; Local Variables: ***
+ ;;; mode: C ***
+ ;;; End: ***
+ */
+
Index: liboctave/Sparse.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/Sparse.cc,v
retrieving revision 1.10
diff -c -r1.10 Sparse.cc
*** liboctave/Sparse.cc 16 Mar 2006 17:48:56 -0000 1.10
--- liboctave/Sparse.cc 2 Apr 2006 20:23:26 -0000
***************
*** 851,857 ****
if (r == dim1 () && c == dim2 ())
return;
! typename Sparse<T>::SparseRep *old_rep = Sparse<T>::rep;
octave_idx_type nc = cols ();
octave_idx_type nr = rows ();
--- 851,858 ----
if (r == dim1 () && c == dim2 ())
return;
! typename Sparse<T>::SparseRep *old_rep = rep;
!
octave_idx_type nc = cols ();
octave_idx_type nr = rows ();
***************
*** 866,900 ****
if (r >= nr)
{
if (c > nc)
! n = cidx(nc);
else
! n = cidx(c);
tmpval = Sparse<T> (r, c, n);
if (c > nc)
{
for (octave_idx_type i = 0; i < nc; i++)
! tmpval.cidx(i) = cidx(i);
for (octave_idx_type i = nc+2; i < c; i++)
tmpval.cidx(i) = tmpval.cidx(i-1);
}
else if (c <= nc)
for (octave_idx_type i = 0; i < c; i++)
! tmpval.cidx(i) = cidx(i);
for (octave_idx_type i = 0; i < n; i++)
{
! tmpval.data(i) = data(i);
! tmpval.ridx(i) = ridx(i);
}
}
else
{
// Count how many non zero terms before we do anything
for (octave_idx_type i = 0; i < c; i++)
! for (octave_idx_type j = cidx(i); j < cidx(i+1); j++)
! if (ridx(j) < r)
n++;
if (n)
--- 867,901 ----
if (r >= nr)
{
if (c > nc)
! n = xcidx(nc);
else
! n = xcidx(c);
tmpval = Sparse<T> (r, c, n);
if (c > nc)
{
for (octave_idx_type i = 0; i < nc; i++)
! tmpval.cidx(i) = xcidx(i);
for (octave_idx_type i = nc+2; i < c; i++)
tmpval.cidx(i) = tmpval.cidx(i-1);
}
else if (c <= nc)
for (octave_idx_type i = 0; i < c; i++)
! tmpval.cidx(i) = xcidx(i);
for (octave_idx_type i = 0; i < n; i++)
{
! tmpval.data(i) = xdata(i);
! tmpval.ridx(i) = xridx(i);
}
}
else
{
// Count how many non zero terms before we do anything
for (octave_idx_type i = 0; i < c; i++)
! for (octave_idx_type j = xcidx(i); j < xcidx(i+1); j++)
! if (xridx(j) < r)
n++;
if (n)
***************
*** 905,915 ****
tmpval.cidx(0);
for (octave_idx_type i = 0, ii = 0; i < c; i++)
{
! for (octave_idx_type j = cidx(i); j < cidx(i+1); j++)
! if (ridx(j) < r)
{
! tmpval.data(ii) = data(j);
! tmpval.ridx(ii++) = ridx(j);
}
tmpval.cidx(i+1) = ii;
}
--- 906,916 ----
tmpval.cidx(0);
for (octave_idx_type i = 0, ii = 0; i < c; i++)
{
! for (octave_idx_type j = xcidx(i); j < xcidx(i+1); j++)
! if (xridx(j) < r)
{
! tmpval.data(ii) = xdata(j);
! tmpval.ridx(ii++) = xridx(j);
}
tmpval.cidx(i+1) = ii;
}
Index: scripts/general/tril.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/general/tril.m,v
retrieving revision 1.26
diff -c -r1.26 tril.m
*** scripts/general/tril.m 15 Mar 2006 03:11:35 -0000 1.26
--- scripts/general/tril.m 2 Apr 2006 20:23:27 -0000
***************
*** 68,78 ****
if (nargin > 0)
[nr, nc] = size (x);
! if (isa (x, "sparse"))
! retval = sparse (nr, nc);
! else
! retval = zeros (nr, nc, class (x));
! endif
endif
if (nargin == 1)
--- 68,74 ----
if (nargin > 0)
[nr, nc] = size (x);
! retval = resize (resize (x, 0), nr, nc);
endif
if (nargin == 1)
Index: scripts/general/triu.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/general/triu.m,v
retrieving revision 1.20
diff -c -r1.20 triu.m
*** scripts/general/triu.m 9 Feb 2006 09:12:03 -0000 1.20
--- scripts/general/triu.m 2 Apr 2006 20:23:27 -0000
***************
*** 28,38 ****
if (nargin > 0)
[nr, nc] = size (x);
! if (isa (x, "sparse"))
! retval = sparse (nr, nc);
! else
! retval = zeros (nr, nc, class (x));
! endif
endif
if (nargin == 1)
--- 28,34 ----
if (nargin > 0)
[nr, nc] = size (x);
! retval = resize (resize (x, 0), nr, nc);
endif
if (nargin == 1)
Index: scripts/special-matrix/hankel.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/special-matrix/hankel.m,v
retrieving revision 1.24
diff -c -r1.24 hankel.m
*** scripts/special-matrix/hankel.m 6 Mar 2006 21:26:50 -0000 1.24
--- scripts/special-matrix/hankel.m 2 Apr 2006 20:23:27 -0000
***************
*** 51,57 ****
function retval = hankel (c, r)
if (nargin == 1)
! r = zeros (size (c));
elseif (nargin != 2)
usage ("hankel (c, r)");
endif
--- 51,57 ----
function retval = hankel (c, r)
if (nargin == 1)
! r = resize (resize (c, 0), size(c));
elseif (nargin != 2)
usage ("hankel (c, r)");
endif
***************
*** 84,90 ****
## This should probably be done with the colon operator...
! retval = zeros (nr, nc);
for i = 1:min (nr, nc)
retval (1:nr-i+1, i) = c (i:nr);
--- 84,90 ----
## This should probably be done with the colon operator...
! retval = resize (resize (c, 0), nr, nc);
for i = 1:min (nr, nc)
retval (1:nr-i+1, i) = c (i:nr);
***************
*** 100,102 ****
--- 100,108 ----
endfor
endfunction
+
+ %!assert(hankel(1:3),[1,2,3;2,3,0;3,0,0])
+ %!assert(hankel(1),[1]);
+ %!assert(hankel(1:3,3:6),[1,2,3,4;2,3,4,5;3,4,5,6]);
+ %!assert(hankel(1:3,3:4),[1,2;2,3;3,4]);
+ %!assert(hankel(1:3,4:6),[1,2,3;2,3,5;3,5,6]);
Index: scripts/special-matrix/toeplitz.m
===================================================================
RCS file: /usr/local/cvsroot/octave/scripts/special-matrix/toeplitz.m,v
retrieving revision 1.23
diff -c -r1.23 toeplitz.m
*** scripts/special-matrix/toeplitz.m 6 Mar 2006 21:26:50 -0000 1.23
--- scripts/special-matrix/toeplitz.m 2 Apr 2006 20:23:27 -0000
***************
*** 97,103 ****
nc = length (r);
nr = length (c);
! retval = zeros (nr, nc);
for i = 1:min (nc, nr)
retval (i:nr, i) = c (1:nr-i+1);
--- 97,103 ----
nc = length (r);
nr = length (c);
! retval = resize (resize (c, 0), nr, nc);
for i = 1:min (nc, nr)
retval (i:nr, i) = c (1:nr-i+1);
Index: src/data.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/data.cc,v
retrieving revision 1.153
diff -c -r1.153 data.cc
*** src/data.cc 30 Mar 2006 21:34:53 -0000 1.153
--- src/data.cc 2 Apr 2006 20:23:28 -0000
***************
*** 1863,1868 ****
--- 1863,1919 ----
return retval;
}
+ DEFUN (resize, args, ,
+ "-*- texinfo -*-\n\
+ @deftypefn {Built-in Function} {} resize (@var{x}, @var{m})\n\
+ @deftypefnx {Built-in Function} {} resize (@var{x}, @var{m}, @var{n})\n\
+ Resize @var{x} to be dimension @address@hidden where @var{m}\n\
+ and @var{n} are scalar. If @var{n} is not supplied, then resize @var{x}\n\
+ to be dimension @address@hidden if @var{m} is scalar, or if\n\
+ @var{m} is a vector, then the values of this vector respresent the\n\
+ dimensions of the resized matrix.\n\
+ @end deftypefn")
+ {
+ octave_value retval;
+ int nargin = args.length ();
+
+ if (nargin == 2)
+ {
+ Array<double> vec = args(1).vector_value ();
+ int ndim = vec.length ();
+ if (ndim == 1)
+ {
+ octave_idx_type m = static_cast<octave_idx_type> (vec(0));
+ retval = args(0);
+ retval = retval.resize (dim_vector (m, m), true);
+ }
+ else
+ {
+ dim_vector dv;
+ dv.resize (ndim);
+ for (int i = 0; i < ndim; i++)
+ dv(i) = static_cast<octave_idx_type> (vec(i));
+ retval = args(0);
+ retval = retval.resize (dv, true);
+ }
+ }
+ else if (nargin == 3)
+ {
+ octave_idx_type m = static_cast<octave_idx_type>
+ (args(1).scalar_value());
+ octave_idx_type n = static_cast<octave_idx_type>
+ (args(2).scalar_value());
+ if (!error_state)
+ {
+ retval = args(0);
+ retval = retval.resize (dim_vector (m, n), true);
+ }
+ }
+ else
+ print_usage ("resize");
+ return retval;
+ }
+
DEFUN (reshape, args, ,
"-*- texinfo -*-\n\
@deftypefn {Function File} {} reshape (@var{a}, @var{m}, @var{n}, @dots{})\n\
Index: src/oct-map.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/oct-map.cc,v
retrieving revision 1.48
diff -c -r1.48 oct-map.cc
*** src/oct-map.cc 12 Jan 2006 20:57:27 -0000 1.48
--- src/oct-map.cc 2 Apr 2006 20:23:28 -0000
***************
*** 141,147 ****
}
Octave_map
! Octave_map::resize (const dim_vector& dv) const
{
Octave_map retval;
--- 141,147 ----
}
Octave_map
! Octave_map::resize (const dim_vector& dv, bool fill) const
{
Octave_map retval;
***************
*** 150,156 ****
for (const_iterator p = begin (); p != end (); p++)
{
Cell tmp = contents(p);
! tmp.resize(dv);
retval.assign (key(p), tmp);
}
--- 150,159 ----
for (const_iterator p = begin (); p != end (); p++)
{
Cell tmp = contents(p);
! if (fill)
! tmp.resize(dv, Cell::resize_fill_value ());
! else
! tmp.resize(dv);
retval.assign (key(p), tmp);
}
Index: src/oct-map.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/oct-map.h,v
retrieving revision 1.43
diff -c -r1.43 oct-map.h
*** src/oct-map.h 12 Jan 2006 17:55:22 -0000 1.43
--- src/oct-map.h 2 Apr 2006 20:23:28 -0000
***************
*** 127,133 ****
Octave_map reshape (const dim_vector& new_dims) const;
! Octave_map resize (const dim_vector& dv) const;
octave_idx_type numel (void) const;
--- 127,133 ----
Octave_map reshape (const dim_vector& new_dims) const;
! Octave_map resize (const dim_vector& dv, bool fill = false) const;
octave_idx_type numel (void) const;
Index: src/ov-base-mat.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-base-mat.cc,v
retrieving revision 1.35
diff -c -r1.35 ov-base-mat.cc
*** src/ov-base-mat.cc 16 Nov 2005 18:45:32 -0000 1.35
--- src/ov-base-mat.cc 2 Apr 2006 20:23:28 -0000
***************
*** 202,207 ****
--- 202,219 ----
}
template <class MT>
+ octave_value
+ octave_base_matrix<MT>::resize (const dim_vector& dv, bool fill) const
+ {
+ MT retval (matrix);
+ if (fill)
+ retval.resize (dv, 0);
+ else
+ retval.resize (dv);
+ return retval;
+ }
+
+ template <class MT>
bool
octave_base_matrix<MT>::is_true (void) const
{
Index: src/ov-base-mat.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-base-mat.h,v
retrieving revision 1.35
diff -c -r1.35 ov-base-mat.h
*** src/ov-base-mat.h 31 Jan 2006 03:43:41 -0000 1.35
--- src/ov-base-mat.h 2 Apr 2006 20:23:28 -0000
***************
*** 102,109 ****
octave_value permute (const Array<int>& vec, bool inv = false) const
{ return MT (matrix.permute (vec, inv)); }
! octave_value resize (const dim_vector& dv) const
! { MT retval (matrix); retval.resize (dv); return retval; }
octave_value all (int dim = 0) const { return matrix.all (dim); }
octave_value any (int dim = 0) const { return matrix.any (dim); }
--- 102,108 ----
octave_value permute (const Array<int>& vec, bool inv = false) const
{ return MT (matrix.permute (vec, inv)); }
! octave_value resize (const dim_vector& dv, bool fill = false) const;
octave_value all (int dim = 0) const { return matrix.all (dim); }
octave_value any (int dim = 0) const { return matrix.any (dim); }
Index: src/ov-base-sparse.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-base-sparse.cc,v
retrieving revision 1.8
diff -c -r1.8 ov-base-sparse.cc
*** src/ov-base-sparse.cc 31 Jan 2006 18:23:00 -0000 1.8
--- src/ov-base-sparse.cc 2 Apr 2006 20:23:28 -0000
***************
*** 40,45 ****
--- 40,46 ----
#include "boolSparse.h"
#include "ov-base-sparse.h"
+ #include "pager.h"
template <class T>
octave_value
***************
*** 193,198 ****
--- 194,207 ----
typ.invalidate_type ();
}
+ template <class T>
+ octave_value
+ octave_base_sparse<T>::resize (const dim_vector& dv, bool) const
+ {
+ T retval (matrix);
+ retval.resize (dv);
+ return retval;
+ }
template <class T>
bool
Index: src/ov-base-sparse.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-base-sparse.h,v
retrieving revision 1.7
diff -c -r1.7 ov-base-sparse.h
*** src/ov-base-sparse.h 22 Feb 2006 20:15:07 -0000 1.7
--- src/ov-base-sparse.h 2 Apr 2006 20:23:28 -0000
***************
*** 112,119 ****
octave_value permute (const Array<int>& vec, bool inv = false) const
{ return T (matrix.permute (vec, inv)); }
! octave_value resize (const dim_vector& dv) const
! { T retval (matrix); retval.resize (dv); return retval; }
octave_value all (int dim = 0) const { return matrix.all (dim); }
octave_value any (int dim = 0) const { return matrix.any (dim); }
--- 112,118 ----
octave_value permute (const Array<int>& vec, bool inv = false) const
{ return T (matrix.permute (vec, inv)); }
! octave_value resize (const dim_vector& dv, bool = false) const;
octave_value all (int dim = 0) const { return matrix.all (dim); }
octave_value any (int dim = 0) const { return matrix.any (dim); }
Index: src/ov-base.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-base.cc,v
retrieving revision 1.75
diff -c -r1.75 ov-base.cc
*** src/ov-base.cc 24 Mar 2006 16:42:44 -0000 1.75
--- src/ov-base.cc 2 Apr 2006 20:23:29 -0000
***************
*** 202,208 ****
}
octave_value
! octave_base_value::resize (const dim_vector&) const
{
gripe_wrong_type_arg ("octave_base_value::resize ()", type_name ());
return octave_value ();
--- 202,208 ----
}
octave_value
! octave_base_value::resize (const dim_vector&, bool) const
{
gripe_wrong_type_arg ("octave_base_value::resize ()", type_name ());
return octave_value ();
Index: src/ov-base.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-base.h,v
retrieving revision 1.86
diff -c -r1.86 ov-base.h
*** src/ov-base.h 24 Mar 2006 16:42:44 -0000 1.86
--- src/ov-base.h 2 Apr 2006 20:23:29 -0000
***************
*** 107,113 ****
octave_value permute (const Array<int>& vec, bool = false) const;
! octave_value resize (const dim_vector&) const;
bool is_defined (void) const { return false; }
--- 107,113 ----
octave_value permute (const Array<int>& vec, bool = false) const;
! octave_value resize (const dim_vector&, bool fill = false) const;
bool is_defined (void) const { return false; }
Index: src/ov-bool.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-bool.cc,v
retrieving revision 1.24
diff -c -r1.24 ov-bool.cc
*** src/ov-bool.cc 26 Apr 2005 19:24:33 -0000 1.24
--- src/ov-bool.cc 2 Apr 2006 20:23:29 -0000
***************
*** 103,108 ****
--- 103,127 ----
return retval;
}
+ octave_value
+ octave_bool::resize (const dim_vector& dv, bool fill) const
+ {
+ if (fill)
+ {
+ boolNDArray retval (dv, false);
+ if (dv.numel())
+ retval(0) = scalar;
+ return retval;
+ }
+ else
+ {
+ boolNDArray retval (dv);
+ if (dv.numel())
+ retval(0) = scalar;
+ return retval;
+ }
+ }
+
octave_value
octave_bool::convert_to_str_internal (bool, bool, char type) const
{
Index: src/ov-bool.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-bool.h,v
retrieving revision 1.30
diff -c -r1.30 ov-bool.h
*** src/ov-bool.h 10 Nov 2005 21:40:48 -0000 1.30
--- src/ov-bool.h 2 Apr 2006 20:23:29 -0000
***************
*** 149,156 ****
boolNDArray bool_array_value (void) const
{ return boolNDArray (dim_vector (1, 1), scalar); }
! octave_value resize (const dim_vector& dv) const
! { boolNDArray retval (dv); if (dv.numel()) retval(0) = scalar; return
retval; }
octave_value convert_to_str_internal (bool pad, bool force, char type)
const;
--- 149,155 ----
boolNDArray bool_array_value (void) const
{ return boolNDArray (dim_vector (1, 1), scalar); }
! octave_value resize (const dim_vector& dv, bool fill = false) const;
octave_value convert_to_str_internal (bool pad, bool force, char type)
const;
Index: src/ov-complex.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-complex.cc,v
retrieving revision 1.30
diff -c -r1.30 ov-complex.cc
*** src/ov-complex.cc 15 Sep 2005 15:36:26 -0000 1.30
--- src/ov-complex.cc 2 Apr 2006 20:23:29 -0000
***************
*** 149,154 ****
--- 149,177 ----
return ComplexNDArray (dim_vector (1, 1), scalar);
}
+ octave_value
+ octave_complex::resize (const dim_vector& dv, bool fill) const
+ {
+ if (fill)
+ {
+ ComplexNDArray retval (dv, ComplexNDArray::resize_fill_value ());
+
+ if (dv.numel ())
+ retval(0) = scalar;
+
+ return retval;
+ }
+ else
+ {
+ ComplexNDArray retval (dv);
+
+ if (dv.numel ())
+ retval(0) = scalar;
+
+ return retval;
+ }
+ }
+
bool
octave_complex::save_ascii (std::ostream& os, bool& infnan_warned,
bool strip_nan_and_inf)
Index: src/ov-complex.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-complex.h,v
retrieving revision 1.31
diff -c -r1.31 ov-complex.h
*** src/ov-complex.h 26 Apr 2005 19:24:33 -0000 1.31
--- src/ov-complex.h 2 Apr 2006 20:23:29 -0000
***************
*** 101,115 ****
NDArray array_value (bool = false) const;
! octave_value resize (const dim_vector& dv) const
! {
! ComplexNDArray retval (dv);
!
! if (dv.numel ())
! retval(0) = scalar;
!
! return retval;
! }
Complex complex_value (bool = false) const;
--- 101,107 ----
NDArray array_value (bool = false) const;
! octave_value resize (const dim_vector& dv, bool fill = false) const;
Complex complex_value (bool = false) const;
Index: src/ov-intx.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-intx.h,v
retrieving revision 1.13
diff -c -r1.13 ov-intx.h
*** src/ov-intx.h 10 Nov 2005 21:40:48 -0000 1.13
--- src/ov-intx.h 2 Apr 2006 20:23:29 -0000
***************
*** 255,266 ****
uint64_array_value (void) const
{ return uint64NDArray (dim_vector (1, 1), scalar); }
! octave_value resize (const dim_vector& dv) const
{
! OCTAVE_INT_NDARRAY_T retval (dv);
! if (dv.numel())
! retval(0) = scalar;
! return retval;
}
double double_value (bool = false) const { return double (scalar); }
--- 255,276 ----
uint64_array_value (void) const
{ return uint64NDArray (dim_vector (1, 1), scalar); }
! octave_value resize (const dim_vector& dv, bool fill = false) const
{
! if (fill)
! {
! OCTAVE_INT_NDARRAY_T retval (dv, 0);
! if (dv.numel())
! retval(0) = scalar;
! return retval;
! }
! else
! {
! OCTAVE_INT_NDARRAY_T retval (dv);
! if (dv.numel())
! retval(0) = scalar;
! return retval;
! }
}
double double_value (bool = false) const { return double (scalar); }
Index: src/ov-range.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-range.cc,v
retrieving revision 1.43
diff -c -r1.43 ov-range.cc
*** src/ov-range.cc 26 Apr 2005 19:24:33 -0000 1.43
--- src/ov-range.cc 2 Apr 2006 20:23:29 -0000
***************
*** 208,213 ****
--- 208,224 ----
return retval;
}
+ octave_value
+ octave_range::resize (const dim_vector& dv, bool fill) const
+ {
+ NDArray retval = array_value ();
+ if (fill)
+ retval.resize (dv, NDArray::resize_fill_value());
+ else
+ retval.resize (dv);
+ return retval;
+ }
+
octave_value
octave_range::convert_to_str_internal (bool pad, bool force, char type) const
{
Index: src/ov-range.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-range.h,v
retrieving revision 1.49
diff -c -r1.49 ov-range.h
*** src/ov-range.h 30 Aug 2005 19:29:50 -0000 1.49
--- src/ov-range.h 2 Apr 2006 20:23:30 -0000
***************
*** 106,113 ****
return dim_vector (n > 0, n);
}
! octave_value resize (const dim_vector& dv) const
! { NDArray retval = array_value (); retval.resize (dv); return retval; }
size_t byte_size (void) const { return 3 * sizeof (double); }
--- 106,113 ----
return dim_vector (n > 0, n);
}
! octave_value resize (const dim_vector& dv, bool fill = false) const;
!
size_t byte_size (void) const { return 3 * sizeof (double); }
Index: src/ov-scalar.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-scalar.cc,v
retrieving revision 1.31
diff -c -r1.31 ov-scalar.cc
*** src/ov-scalar.cc 26 Apr 2005 19:24:33 -0000 1.31
--- src/ov-scalar.cc 2 Apr 2006 20:23:30 -0000
***************
*** 106,111 ****
--- 106,134 ----
return retval;
}
+ octave_value
+ octave_scalar::resize (const dim_vector& dv, bool fill) const
+ {
+ if (fill)
+ {
+ NDArray retval (dv, NDArray::resize_fill_value());
+
+ if (dv.numel ())
+ retval(0) = scalar;
+
+ return retval;
+ }
+ else
+ {
+ NDArray retval (dv);
+
+ if (dv.numel ())
+ retval(0) = scalar;
+
+ return retval;
+ }
+ }
+
octave_value
octave_scalar::convert_to_str_internal (bool, bool, char type) const
{
Index: src/ov-scalar.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-scalar.h,v
retrieving revision 1.41
diff -c -r1.41 ov-scalar.h
*** src/ov-scalar.h 10 Nov 2005 21:40:48 -0000 1.41
--- src/ov-scalar.h 2 Apr 2006 20:23:30 -0000
***************
*** 138,152 ****
NDArray array_value (bool = false) const
{ return NDArray (dim_vector (1, 1), scalar); }
! octave_value resize (const dim_vector& dv) const
! {
! NDArray retval (dv);
!
! if (dv.numel ())
! retval(0) = scalar;
!
! return retval;
! }
Complex complex_value (bool = false) const { return scalar; }
--- 138,144 ----
NDArray array_value (bool = false) const
{ return NDArray (dim_vector (1, 1), scalar); }
! octave_value resize (const dim_vector& dv, bool fill = false) const;
Complex complex_value (bool = false) const { return scalar; }
Index: src/ov-str-mat.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-str-mat.cc,v
retrieving revision 1.55
diff -c -r1.55 ov-str-mat.cc
*** src/ov-str-mat.cc 24 Mar 2006 16:42:44 -0000 1.55
--- src/ov-str-mat.cc 2 Apr 2006 20:23:30 -0000
***************
*** 149,154 ****
--- 149,165 ----
::assign (matrix, tmp, Vstring_fill_char);
}
+ octave_value
+ octave_char_matrix_str::resize (const dim_vector& dv, bool fill) const
+ {
+ charNDArray retval (matrix);
+ if (fill)
+ retval.resize (dv, charNDArray::resize_fill_value());
+ else
+ retval.resize (dv);
+ return octave_value (retval, true);
+ }
+
bool
octave_char_matrix_str::valid_as_scalar_index (void) const
{
Index: src/ov-str-mat.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-str-mat.h,v
retrieving revision 1.41
diff -c -r1.41 ov-str-mat.h
*** src/ov-str-mat.h 24 Mar 2006 16:42:44 -0000 1.41
--- src/ov-str-mat.h 2 Apr 2006 20:23:30 -0000
***************
*** 96,107 ****
octave_value permute (const Array<int>& vec, bool inv = false) const
{ return octave_value (charNDArray (matrix.permute (vec, inv)), true); }
! octave_value resize (const dim_vector& dv) const
! {
! charNDArray retval (matrix);
! retval.resize (dv);
! return octave_value (retval, true);
! }
bool is_string (void) const { return true; }
--- 96,102 ----
octave_value permute (const Array<int>& vec, bool inv = false) const
{ return octave_value (charNDArray (matrix.permute (vec, inv)), true); }
! octave_value resize (const dim_vector& dv, bool fill = false) const;
bool is_string (void) const { return true; }
Index: src/ov.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov.h,v
retrieving revision 1.132
diff -c -r1.132 ov.h
*** src/ov.h 24 Mar 2006 16:42:44 -0000 1.132
--- src/ov.h 2 Apr 2006 20:23:30 -0000
***************
*** 391,398 ****
octave_value ipermute (const Array<int>& vec) const
{ return rep->permute (vec, true); }
! virtual octave_value resize (const dim_vector& dv) const
! { return rep->resize (dv);}
// Does this constant have a type? Both of these are provided since
// it is sometimes more natural to write is_undefined() instead of
--- 391,398 ----
octave_value ipermute (const Array<int>& vec) const
{ return rep->permute (vec, true); }
! virtual octave_value resize (const dim_vector& dv, bool fill = false) const
! { return rep->resize (dv, fill);}
// Does this constant have a type? Both of these are provided since
// it is sometimes more natural to write is_undefined() instead of
Index: src/DLD-FUNCTIONS/fftw_wisdom.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/DLD-FUNCTIONS/fftw_wisdom.cc,v
retrieving revision 1.7
diff -c -r1.7 fftw_wisdom.cc
*** src/DLD-FUNCTIONS/fftw_wisdom.cc 6 Mar 2006 21:26:54 -0000 1.7
--- src/DLD-FUNCTIONS/fftw_wisdom.cc 1 Apr 2006 13:55:14 -0000
***************
*** 108,125 ****
overwrite = true;
}
std::string wisdom = octave_env::make_absolute
! (Vload_path_dir_path.find_first_of (args(0).string_value ()),
! octave_env::getcwd ());
// XXX FIXME XXX -- should probably protect FILE* resources with
// auto_ptr or similar...
if (wisdom.empty () || overwrite)
{
! FILE *ofile = fopen (wisdom.c_str (), "wb");
! fftw_export_wisdom_to_file (ofile);
! fclose (ofile);
}
else
{
--- 108,135 ----
overwrite = true;
}
+ std::string str = args(0).string_value ();
std::string wisdom = octave_env::make_absolute
! (Vload_path_dir_path.find_first_of (str), octave_env::getcwd ());
// XXX FIXME XXX -- should probably protect FILE* resources with
// auto_ptr or similar...
if (wisdom.empty () || overwrite)
{
! if (str.empty ())
! error ("fftw_wisdom: can not save to file");
! else
! {
! FILE *ofile = fopen (str.c_str (), "wb");
! if (! ofile)
! error ("fftw_wisdom: can not save to file %s", str.c_str());
! else
! {
! fftw_export_wisdom_to_file (ofile);
! fclose (ofile);
! }
! }
}
else
{
- Importing FIXES functions from octave-forge, and a couple of other things,
David Bateman <=
- Re: Importing FIXES functions from octave-forge, and a couple of other things, Paul Kienzle, 2006/04/02
- Importing FIXES functions from octave-forge, and a couple of other things, John W. Eaton, 2006/04/03
- Importing FIXES functions from octave-forge, and a couple of other things, John W. Eaton, 2006/04/03
- Importing FIXES functions from octave-forge, and a couple of other things, John W. Eaton, 2006/04/03
- Importing FIXES functions from octave-forge, and a couple of other things, John W. Eaton, 2006/04/03
- Importing FIXES functions from octave-forge, and a couple of other things, John W. Eaton, 2006/04/03