Mercurial > hg > octave-lyh
changeset 14817:67897baaa05f
Adapt implementation of betaincinv to Octave.
Add support for 'single' type variables.
Use superclass Array rather than Matrix or NDArray in function prototypes.
* lo-specfun.h, lo-specfun.cc: Use superclass Array rather than Matrix or NDArray
in function prototypes for betaincinv.
* beta.m: Add seealso links in docstring.
* betainc.cc (betaincinv): Move code from mappers.cc. Cast output
to float if any of the inputs to function are of type single.
* mappers.cc (betaincinv): Delete code for betaincinv.
| author | Rik <octave@nomad.inbox5.com> |
|---|---|
| date | Thu, 28 Jun 2012 10:18:38 -0700 |
| parents | 0a868d90436b |
| children | c6ae30f73946 |
| files | liboctave/lo-specfun.cc liboctave/lo-specfun.h scripts/specfun/beta.m src/DLD-FUNCTIONS/betainc.cc src/mappers.cc |
| diffstat | 5 files changed, 269 insertions(+), 364 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/lo-specfun.cc +++ b/liboctave/lo-specfun.cc @@ -3153,7 +3153,7 @@ // Algorithm based on the one by John Burkardt. // See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html // -// The original code is distributed under the GNU LGPL license. +// The original code is distributed under the GNU LGPL v3 license. // // Reference: // @@ -3269,7 +3269,7 @@ // Algorithm based on the one by John Burkardt. // See http://people.sc.fsu.edu/~jburkardt/cpp_src/asa109/asa109.html // -// The original code is distributed under the GNU LGPL license. +// The original code is distributed under the GNU LGPL v3 license. // // Reference: // @@ -3472,150 +3472,142 @@ return value; } -Matrix -betaincinv (double x, double a, const Matrix& b) +Array<double> +betaincinv (double x, double a, const Array<double>& b) { - octave_idx_type nr = b.rows (); - octave_idx_type nc = b.cols (); - - Matrix retval (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betaincinv (x, a, b(i,j)); + dim_vector dv = b.dims (); + octave_idx_type nel = dv.numel (); + + Array<double> retval (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x, a, b(i)); return retval; } -Matrix -betaincinv (double x, const Matrix& a, double b) +Array<double> +betaincinv (double x, const Array<double>& a, double b) { - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - - Matrix retval (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betaincinv (x, a(i,j), b); + dim_vector dv = a.dims (); + octave_idx_type nel = dv.numel (); + + Array<double> retval (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x, a(i), b); return retval; } -Matrix -betaincinv (double x, const Matrix& a, const Matrix& b) +Array<double> +betaincinv (double x, const Array<double>& a, const Array<double>& b) { - Matrix retval; - - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - - if (a_nr == b_nr && a_nc == b_nc) + Array<double> retval; + dim_vector dv = a.dims (); + + if (dv == b.dims ()) { - retval.resize (a_nr, a_nc); - - for (octave_idx_type j = 0; j < a_nc; j++) - for (octave_idx_type i = 0; i < a_nr; i++) - retval(i,j) = betaincinv (x, a(i,j), b(i,j)); + octave_idx_type nel = dv.numel (); + + retval.resize (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x, a(i), b(i)); } else - gripe_betaincinv_nonconformant (1, 1, a_nr, a_nc, b_nr, b_nc); + gripe_betaincinv_nonconformant (dim_vector (0, 0), dv, b.dims ()); return retval; } -Matrix -betaincinv (const Matrix& x, double a, double b) +Array<double> +betaincinv (const Array<double>& x, double a, double b) { - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - Matrix retval (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betaincinv (x(i,j), a, b); + dim_vector dv = x.dims (); + octave_idx_type nel = dv.numel (); + + Array<double> retval (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x(i), a, b); return retval; } -Matrix -betaincinv (const Matrix& x, double a, const Matrix& b) +Array<double> +betaincinv (const Array<double>& x, double a, const Array<double>& b) { - Matrix retval; - - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - - if (nr == b_nr && nc == b_nc) + Array<double> retval; + dim_vector dv = x.dims (); + + if (dv == b.dims ()) { - retval.resize (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betaincinv (x(i,j), a, b(i,j)); + octave_idx_type nel = dv.numel (); + + retval.resize (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x(i), a, b(i)); } else - gripe_betaincinv_nonconformant (nr, nc, 1, 1, b_nr, b_nc); + gripe_betaincinv_nonconformant (dv, dim_vector (0, 0), b.dims ()); return retval; } -Matrix -betaincinv (const Matrix& x, const Matrix& a, double b) +Array<double> +betaincinv (const Array<double>& x, const Array<double>& a, double b) { - Matrix retval; - - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - if (nr == a_nr && nc == a_nc) + Array<double> retval; + dim_vector dv = x.dims (); + + if (dv == a.dims ()) { - retval.resize (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betaincinv (x(i,j), a(i,j), b); + octave_idx_type nel = dv.numel (); + + retval.resize (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x(i), a(i), b); } else - gripe_betaincinv_nonconformant (nr, nc, a_nr, a_nc, 1, 1); + gripe_betaincinv_nonconformant (dv, a.dims (), dim_vector (0, 0)); return retval; } -Matrix -betaincinv (const Matrix& x, const Matrix& a, const Matrix& b) +Array<double> +betaincinv (const Array<double>& x, const Array<double>& a, const Array<double>& b) { - Matrix retval; - - octave_idx_type nr = x.rows (); - octave_idx_type nc = x.cols (); - - octave_idx_type a_nr = a.rows (); - octave_idx_type a_nc = a.cols (); - - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - - if (nr == a_nr && nr == b_nr && nc == a_nc && nc == b_nc) + Array<double> retval; + dim_vector dv = x.dims (); + + if (dv == a.dims () && dv == b.dims ()) { - retval.resize (nr, nc); - - for (octave_idx_type j = 0; j < nc; j++) - for (octave_idx_type i = 0; i < nr; i++) - retval(i,j) = betaincinv (x(i,j), a(i,j), b(i,j)); + octave_idx_type nel = dv.numel (); + + retval.resize (dv); + + double *pretval = retval.fortran_vec (); + + for (octave_idx_type i = 0; i < nel; i++) + *pretval++ = betaincinv (x(i), a(i), b(i)); } else - gripe_betaincinv_nonconformant (nr, nc, a_nr, a_nc, b_nr, b_nc); + gripe_betaincinv_nonconformant (dv, a.dims (), b.dims ()); return retval; }
--- a/liboctave/lo-specfun.h +++ b/liboctave/lo-specfun.h @@ -581,16 +581,14 @@ extern OCTAVE_API double erfcx (double x); extern OCTAVE_API float erfcx (float x); -extern OCTAVE_API double betaincinv (double y, double p, double q); -extern OCTAVE_API Matrix betaincinv (double x, double a, const Matrix& b); -extern OCTAVE_API Matrix betaincinv (double x, const Matrix& a, double b); -extern OCTAVE_API Matrix betaincinv (double x, const Matrix& a, const Matrix& b); - -extern OCTAVE_API Matrix betaincinv (const Matrix& x, double a, double b); -extern OCTAVE_API Matrix betaincinv (const Matrix& x, double a, const Matrix& b); -extern OCTAVE_API Matrix betaincinv (const Matrix& x, const Matrix& a, double b); -extern OCTAVE_API Matrix betaincinv (const Matrix& x, const Matrix& a, const Matrix& b); - -extern OCTAVE_API double betain (double x, double p, double q, double beta, bool& err); +extern OCTAVE_API double betaincinv (double x, double a, double b); +extern OCTAVE_API Array<double> betaincinv (double x, double a, const Array<double>& b); +extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, double b); +extern OCTAVE_API Array<double> betaincinv (double x, const Array<double>& a, const Array<double>& b); +extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, double b); +extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, double b); +extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, double a, const Array<double>& b); +extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, double b); +extern OCTAVE_API Array<double> betaincinv (const Array<double>& x, const Array<double>& a, const Array<double>& b); #endif
--- a/scripts/specfun/beta.m +++ b/scripts/specfun/beta.m @@ -31,6 +31,7 @@ ## @end example ## ## @end ifnottex +## @seealso{betaln, betainc} ## @end deftypefn ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
--- a/src/DLD-FUNCTIONS/betainc.cc +++ b/src/DLD-FUNCTIONS/betainc.cc @@ -32,6 +32,9 @@ #include "oct-obj.h" #include "utils.h" +// FIXME: These functions do not need to be dynamically loaded. They should +// be placed elsewhere in the Octave code hierarchy. + DEFUN_DLD (betainc, args, , "-*- texinfo -*-\n\ @deftypefn {Mapping Function} {} betainc (@var{x}, @var{a}, @var{b})\n\ @@ -324,3 +327,164 @@ %!error betainc (1,2) %!error betainc (1,2,3,4) */ + +DEFUN_DLD (betaincinv, args, , + "-*- texinfo -*-\n\ +@deftypefn {Mapping Function} {} betaincinv (@var{y}, @var{a}, @var{b})\n\ +Compute the inverse of the incomplete Beta function, i.e., @var{x} such that\n\ +\n\ +@example\n\ +@var{y} == betainc (@var{x}, @var{a}, @var{b}) \n\ +@end example\n\ +@seealso{betainc, beta, betaln}\n\ +@end deftypefn") +{ + octave_value retval; + + int nargin = args.length (); + + if (nargin == 3) + { + octave_value x_arg = args(0); + octave_value a_arg = args(1); + octave_value b_arg = args(2); + + if (x_arg.is_scalar_type ()) + { + double x = x_arg.double_value (); + + if (a_arg.is_scalar_type ()) + { + double a = a_arg.double_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + Array<double> b = b_arg.array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + else + { + Array<double> a = a_arg.array_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + Array<double> b = b_arg.array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + } + else + { + Array<double> x = x_arg.array_value (); + + if (a_arg.is_scalar_type ()) + { + double a = a_arg.double_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + Array<double> b = b_arg.array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + else + { + Array<double> a = a_arg.array_value (); + + if (! error_state) + { + if (b_arg.is_scalar_type ()) + { + double b = b_arg.double_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + else + { + Array<double> b = b_arg.array_value (); + + if (! error_state) + retval = betaincinv (x, a, b); + } + } + } + } + + // FIXME: It would be better to have an algorithm for betaincinv which accepted + // float inputs and returned float outputs. As it is, we do extra work + // to calculate betaincinv to double precision and then throw that precision + // away. + if (x_arg.is_single_type () || a_arg.is_single_type () || + b_arg.is_single_type ()) + { + retval = Array<float> (retval.array_value ()); + } + } + else + print_usage (); + + return retval; +} + +/* +%!assert (betaincinv ([0.875 0.6875], [1 2], 3), [0.5 0.5], sqrt (eps)) +%!assert (betaincinv (0.5, 3, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.34375, 4, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.2265625, 5, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.14453125, 6, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.08984375, 7, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.0546875, 8, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.03271484375, 9, 3), 0.5, sqrt (eps)) +%!assert (betaincinv (0.019287109375, 10, 3), 0.5, sqrt (eps)) + +## Test class single as well +%!assert (betaincinv ([0.875 0.6875], [1 2], single (3)), [0.5 0.5], sqrt (eps ("single"))) +%!assert (betaincinv (0.5, 3, single (3)), 0.5, sqrt (eps ("single"))) +%!assert (betaincinv (0.34375, 4, single (3)), 0.5, sqrt (eps ("single"))) + +## Extreme values +%!assert (betaincinv (0, 42, 42), 0, sqrt (eps)) +%!assert (betaincinv (1, 42, 42), 1, sqrt (eps)) + +%!error betaincinv () +%!error betaincinv (1, 2) +*/ +
--- a/src/mappers.cc +++ b/src/mappers.cc @@ -1945,256 +1945,6 @@ %!error tanh (1, 2) */ -DEFUN (betaincinv, args, , - "-*- texinfo -*-\n\ -@deftypefn {Mapping Function} {} betaincinv (@var{y}, @var{z}, @var{w})\n\ -Compute the inverse of the incomplete Beta function, i.e., @var{x} such that\n\ -\n\ -@example\n\ -@var{y} == betainc (@var{x}, @var{z}, @var{w}) \n\ -@end example\n\ -@seealso{beta, betainc, betaln}\n\ -@end deftypefn") -{ - octave_value retval; - - int nargin = args.length (); - - if (nargin == 3) - { - octave_value x_arg = args(0); - octave_value a_arg = args(1); - octave_value b_arg = args(2); - - // FIXME Can we make a template version of the duplicated code below - if (x_arg.is_single_type () || a_arg.is_single_type () || - b_arg.is_single_type ()) - { - if (x_arg.is_scalar_type ()) - { - float x = x_arg.float_value (); - - if (a_arg.is_scalar_type ()) - { - float a = a_arg.float_value (); - - if (! error_state) - { - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - else - { - FloatNDArray b = b_arg.float_array_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - } - } - else - { - FloatNDArray a = a_arg.float_array_value (); - - if (! error_state) - { - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - else - { - FloatNDArray b = b_arg.float_array_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - } - } - } - else - { - FloatNDArray x = x_arg.float_array_value (); - - if (a_arg.is_scalar_type ()) - { - float a = a_arg.float_value (); - - if (! error_state) - { - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - else - { - FloatNDArray b = b_arg.float_array_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - } - } - else - { - FloatNDArray a = a_arg.float_array_value (); - - if (! error_state) - { - if (b_arg.is_scalar_type ()) - { - float b = b_arg.float_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - else - { - FloatNDArray b = b_arg.float_array_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - } - } - } - } - else - { - if (x_arg.is_scalar_type ()) - { - double x = x_arg.double_value (); - - if (a_arg.is_scalar_type ()) - { - double a = a_arg.double_value (); - - if (! error_state) - { - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - else - { - NDArray b = b_arg.array_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - } - } - else - { - NDArray a = a_arg.array_value (); - - if (! error_state) - { - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - else - { - NDArray b = b_arg.array_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - } - } - } - else - { - NDArray x = x_arg.array_value (); - - if (a_arg.is_scalar_type ()) - { - double a = a_arg.double_value (); - - if (! error_state) - { - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - else - { - NDArray b = b_arg.array_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - } - } - else - { - NDArray a = a_arg.array_value (); - - if (! error_state) - { - if (b_arg.is_scalar_type ()) - { - double b = b_arg.double_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - else - { - NDArray b = b_arg.array_value (); - - if (! error_state) - retval = betaincinv (x, a, b); - } - } - } - } - } - } - else - print_usage (); - - return retval; -} - -/* -%!assert (betaincinv ([0.875 0.6875], [1 2], 3), [0.5 0.5], sqrt (eps)); -%!assert (betaincinv (0.5, 3, 3), 0.5, sqrt (eps)); -%!assert (betaincinv (0.34375, 4, 3), 0.5, sqrt (eps)); -%!assert (betaincinv (0.2265625, 5, 3), 0.5, sqrt (eps)); -%!assert (betaincinv (0.14453125, 6, 3), 0.5, sqrt (eps)); -%!assert (betaincinv (0.08984375, 7, 3), 0.5, sqrt (eps)); -%!assert (betaincinv (0.0546875, 8, 3), 0.5, sqrt (eps)); -%!assert (betaincinv (0.03271484375, 9, 3), 0.5, sqrt (eps)); -%!assert (betaincinv (0.019287109375, 10, 3), 0.5, sqrt (eps)); - -%!assert (betaincinv (0, 42, 42), 0, sqrt (eps)); -%!assert (betaincinv (1, 42, 42), 1, sqrt (eps)); - -%!error betaincinv() -%!error betaincinv(1, 2) -*/ - DEFUNX ("toascii", Ftoascii, args, , "-*- texinfo -*-\n\ @deftypefn {Mapping Function} {} toascii (@var{s})\n\