[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Five other functions that are in Matlab core ported to Octave
|
From: |
David Bateman |
|
Subject: |
Five other functions that are in Matlab core ported to Octave |
|
Date: |
Thu, 06 Sep 2007 14:01:14 +0200 |
|
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
The attached patch contains versions of the typecast, swapbytes, gzip,
celldisp and bsxfun functions for Octave. Most of these functions were
easy enough to write though typecast needed to be implemented as an
oct-file and gzip needed to be slightly more complex than it could have
been as the matlab gzip keeps the original file, and I had to deal with
path issues under windows.
However, the bsxfun function was a little more complex. Basically it is
easy to implement the bsxfun using the octave_value::subsassgn method,
however that is extremely slow. Using data.cc(do_cat_op) was 50 times
faster, but still slower than an implementation using repmat. I finally
special cased most types and used the Array<T>::insert method directly.
So bsxfun has a lot of special casing..
This however points to one issue, ie the huge difference in speed
between octave_value::subsassgn and data.cc(do_cat_op). It seems to be
that a case like
a = zeros(m, n)
for i = 1: n;
a(:,i) = ...
end
could be significantly sped up if the octave_value::subsasgn method
could special case to do_cat_op early. It is unclear to me why the
subsassgn method is so slow, and some profiling should be done, but it
might have something to do this the function dereferencing the the
function that finally does the work. The above subsassgn has a call tree
something like
octave_value::subsassgn
-> octave_matrix::subsassgn
-> octave_base_value::numeric_assign
-> op-m-m.cc (assign)
-> octave_matrix::asssign
-> ::assign
-> ::assign2
perhaps that is just too deep and we are paying for it.. When I get some
time I might profile the above case and try and identify where the slow
up is..
Regards
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
*** ./doc/interpreter/system.txi.orig22 2007-09-05 17:57:19.179773220 +0200
--- ./doc/interpreter/system.txi 2007-09-05 17:57:21.513657364 +0200
***************
*** 195,200 ****
--- 195,202 ----
@DOCSTRING(bunzip2)
+ @DOCSTRING(gzip)
+
@DOCSTRING(gunzip)
@DOCSTRING(tar)
*** ./doc/interpreter/container.txi.orig22 2007-09-05 17:57:47.183382985
+0200
--- ./doc/interpreter/container.txi 2007-09-05 22:39:10.624956969 +0200
***************
*** 482,487 ****
--- 482,493 ----
@end group
@end example
+ In general nested cell arrays are displayed hierarchically as above. In
+ some circumstances it makes sense to reference them by their index, and
+ this can be performed by the @code{celldisp} function.
+
+ @DOCSTRING(celldisp)
+
@menu
* Creating Cell Arrays::
* Indexing Cell Arrays::
*** ./doc/interpreter/matrix.txi.orig22 2007-09-05 22:19:42.363448639 +0200
--- ./doc/interpreter/matrix.txi 2007-09-05 22:20:15.199776507 +0200
***************
*** 144,149 ****
--- 144,151 ----
@DOCSTRING(arrayfun)
+ @DOCSTRING(bsxfun)
+
@node Special Utility Matrices
@section Special Utility Matrices
*** ./doc/interpreter/data.txi.orig22 2007-09-05 17:50:04.786306722 +0200
--- ./doc/interpreter/data.txi 2007-09-05 17:50:08.496123084 +0200
***************
*** 47,52 ****
--- 47,56 ----
@DOCSTRING(cast)
+ @DOCSTRING(typecast)
+
+ @DOCSTRING(swapbytes)
+
@menu
* Numeric Objects::
* Missing Data::
*** ./liboctave/Array-util.cc.orig22 2007-09-05 17:38:57.011356392 +0200
--- ./liboctave/Array-util.cc 2007-09-05 17:38:03.998986597 +0200
***************
*** 60,69 ****
ra_idx(start_dimension)++;
int n = ra_idx.length () - 1;
for (int i = start_dimension; i < n; i++)
{
! if (ra_idx(i) < dimensions(i))
break;
else
{
--- 60,70 ----
ra_idx(start_dimension)++;
int n = ra_idx.length () - 1;
+ int nda = dimensions.length ();
for (int i = start_dimension; i < n; i++)
{
! if (ra_idx(i) < (i < nda ? dimensions(i) : 1))
break;
else
{
*** ./scripts/general/Makefile.in.orig22 2007-09-05 17:53:08.603202141
+0200
--- ./scripts/general/Makefile.in 2007-09-05 22:39:11.597907424 +0200
***************
*** 22,29 ****
SOURCES = __isequal__.m __splinen__.m accumarray.m arrayfun.m bicubic.m \
bitcmp.m bitget.m bitset.m blkdiag.m cart2pol.m cart2sph.m cell2mat.m \
! circshift.m common_size.m cplxpair.m cumtrapz.m deal.m del2.m diff.m \
! flipdim.m fliplr.m flipud.m gradient.m ind2sub.m int2str.m interp1.m \
interp2.m interp3.m interpn.m interpft.m is_duplicate_entry.m isa.m \
isdefinite.m isdir.m isequal.m isequalwithequalnans.m isscalar.m \
issquare.m issymmetric.m isvector.m logical.m logspace.m lookup.m mod.m \
--- 22,29 ----
SOURCES = __isequal__.m __splinen__.m accumarray.m arrayfun.m bicubic.m \
bitcmp.m bitget.m bitset.m blkdiag.m cart2pol.m cart2sph.m cell2mat.m \
! celldisp.m circshift.m common_size.m cplxpair.m cumtrapz.m deal.m del2.m \
! diff.m flipdim.m fliplr.m flipud.m gradient.m ind2sub.m int2str.m interp1.m
\
interp2.m interp3.m interpn.m interpft.m is_duplicate_entry.m isa.m \
isdefinite.m isdir.m isequal.m isequalwithequalnans.m isscalar.m \
issquare.m issymmetric.m isvector.m logical.m logspace.m lookup.m mod.m \
*** ./scripts/general/celldisp.m.orig22 2007-09-05 17:52:22.723475647 +0200
--- ./scripts/general/celldisp.m 2007-09-05 17:52:18.438687938 +0200
***************
*** 0 ****
--- 1,63 ----
+ ## Copyright (C) 2007 David Bateman
+ ##
+ ## This file is part of Octave.
+ ##
+ ## Octave is free software; you can redistribute it and/or modify it
+ ## under the terms of the GNU General Public License as published by
+ ## the Free Software Foundation; either version 2, or (at your option)
+ ## any later version.
+ ##
+ ## Octave is distributed in the hope that it will be useful, but
+ ## WITHOUT ANY WARRANTY; without even the implied warranty of
+ ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ ## General Public License for more details.
+ ##
+ ## You should have received a copy of the GNU General Public License
+ ## along with Octave; see the file COPYING. If not, write to the Free
+ ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ ## 02110-1301, USA.
+
+ ## -*- texinfo -*-
+ ## @deftypefn {Function File} {} celldisp (@var{c}, @var{name})
+ ## Recursively display the contents of a cell array. By default the values
+ ## are displayed with the name of the variable @var{c}. However, this name
+ ## can be replaced with the variable @var{name}.
+ ## @seealso{disp}
+ ## @end deftypefn
+
+ ## This is ugly, but seems to be what matlab does..
+
+ function celldisp (c, name)
+ if (nargin < 1 || nargin > 2)
+ print_usage ();
+ endif
+
+ if (! iscell (c))
+ error ("celldisp: argument must be a cell array");
+ endif
+
+ if (nargin == 1)
+ name = inputname (1);
+ endif
+
+ for i = 1: numel (c)
+ if (iscell (c{i}))
+ celldisp (c{i}, sprintf ("%s{%s}", name, indices (size (c), i)));
+ else
+ disp (sprintf ("%s{%s} = \n", name, indices (size (c), i)));
+ disp (c{i});
+ disp ("");
+ endif
+ endfor
+ endfunction
+
+ function s = indices (dv, i)
+ if (sum (dv != 1) > 1)
+ c = cell (size (dv));
+ [c{:}] = ind2sub (dv, i);
+ s = sprintf("%i,", c{:});
+ s(end) = [];
+ else
+ s = sprintf("%i", i);
+ endif
+ endfunction
*** ./scripts/miscellaneous/swapbytes.m.orig22 2007-09-05 17:50:38.838620939
+0200
--- ./scripts/miscellaneous/swapbytes.m 2007-09-05 17:50:33.794870658 +0200
***************
*** 0 ****
--- 1,57 ----
+ ## Copyright (C) 2007 David Bateman
+ ##
+ ## This file is part of Octave.
+ ##
+ ## Octave is free software; you can redistribute it and/or modify it
+ ## under the terms of the GNU General Public License as published by
+ ## the Free Software Foundation; either version 2, or (at your option)
+ ## any later version.
+ ##
+ ## Octave is distributed in the hope that it will be useful, but
+ ## WITHOUT ANY WARRANTY; without even the implied warranty of
+ ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ ## General Public License for more details.
+ ##
+ ## You should have received a copy of the GNU General Public License
+ ## along with Octave; see the file COPYING. If not, write to the Free
+ ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ ## 02110-1301, USA.
+
+ ## -*- texinfo -*-
+ ## @deftypefn {Function File} {} swapbytes (@var{x})
+ ## Swaps the byte order on values, converting from little endian to big
+ ## endian and visa-versa. For example
+ ##
+ ## @example
+ ## @group
+ ## swapbytes (uint16 (1:4))
+ ## @result{} [ 256 512 768 1024]
+ ## @end group
+ ## @end example
+ ##
+ ## @seealso{typecast, cast}
+ ## @end deftypefn
+
+ function y = swapbytes (x)
+ if (nargin != 1)
+ print_usage ();
+ endif
+
+ clx = class (x);
+ if (strcmp (clx, "int8") || strcmp (clx, "uint8") || isempty (x))
+ y = x;
+ else
+ if (strcmp (clx, "int16") || strcmp (clx, "uint16"))
+ nb = 2;
+ elseif (strcmp (clx, "int32") || strcmp (clx, "uint32"))
+ nb = 4;
+ elseif (strcmp (clx, "int64") || strcmp (clx, "uint64") ||
+ strcmp (clx, "double"))
+ nb = 8;
+ else
+ error ("swapbytes: invalid class of object");
+ endif
+ y = reshape (typecast (reshape (typecast (x(:), "uint8"), nb, numel (x))
+ ([nb : -1 : 1], :) (:), clx), size(x));
+ endif
+ endfunction
*** ./scripts/miscellaneous/gzip.m.orig22 2007-09-05 17:51:04.761337357
+0200
--- ./scripts/miscellaneous/gzip.m 2007-09-05 17:50:59.659589992 +0200
***************
*** 0 ****
--- 1,110 ----
+ ## Copyright (C) 2007 David Bateman
+ ##
+ ## This file is part of Octave.
+ ##
+ ## Octave is free software; you can redistribute it and/or modify it
+ ## under the terms of the GNU General Public License as published by
+ ## the Free Software Foundation; either version 2, or (at your option)
+ ## any later version.
+ ##
+ ## Octave is distributed in the hope that it will be useful, but
+ ## WITHOUT ANY WARRANTY; without even the implied warranty of
+ ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ ## General Public License for more details.
+ ##
+ ## You should have received a copy of the GNU General Public License
+ ## along with Octave; see the file COPYING. If not, write to the Free
+ ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ ## 02110-1301, USA.
+
+ ## -*- texinfo -*-
+ ## @deftypefn {Function File} address@hidden =} gzip (@var{files})
+ ## @deftypefnx {Function File} address@hidden =} gzip (@var{files},
@var{outdir})
+ ## Compress the list of files and/or directories specified in @var{files}.
+ ## Each file is compressed separately and a new file with a '.gz' extension
+ ## is create. The original file is not touch. If @var{rootdir} is defined
+ ## the compressed versions of the files are placed in this directory.
+ ## @seealso{gunzip, zip, tar}
+ ## @end deftypefn
+
+ function entries = gzip (files, outdir)
+
+ if (nargin == 1 || nargin == 2)
+
+ if (nargin == 2 && ! exist (outdir, "dir"))
+ error ("gzip: output directory does not exist");
+ endif
+
+ if (ischar (files))
+ files = cellstr (files);
+ endif
+
+ if (nargin == 1)
+ outdir = tmpnam ();
+ mkdir (outdir);
+ endif
+
+ cwd = pwd();
+ unwind_protect
+ if (iscellstr (files))
+ files = glob (files);
+
+ ## Ignore any file with a .gz extension
+ files (cellfun (@(x) strcmp (x(end-2:end), ".gz"), files)) = [];
+
+ copyfile (files, outdir);
+ [d, f] = myfileparts(files);
+ cd (outdir);
+
+ cmd = sprintf ("gzip -r %s", sprintf (" %s", f{:}));
+
+ [status, output] = system (cmd);
+
+ if (status == 0)
+
+ if (nargin == 2)
+ gzfiles = cellfun(@(x) fullfile (outdir, sprintf ("%s.gz", x)), ...
+ f, "UniformOutput", false);
+ else
+ movefile (cellfun(@(x) sprintf ("%s.gz", x), f, ...
+ "UniformOutput", false), cwd);
+ gzfiles = cellfun(@(x) sprintf ("%s.gz", x), ...
+ files, "UniformOutput", false);
+ endif
+
+ if (nargout > 0)
+ entries = gzfiles;
+ endif
+ else
+ error ("gzip: failed with exit status = %d", status);
+ endif
+
+ else
+ error ("gzip: expecting all arguments to be character strings");
+ endif
+ unwind_protect_cleanup
+ cd(cwd);
+ if (nargin == 1)
+ crr = confirm_recursive_rmdir ();
+ unwind_protect
+ confirm_recursive_rmdir (false);
+ rmdir (outdir, "s");
+ unwind_protect_cleanup
+ confirm_recursive_rmdir (crr);
+ end_unwind_protect
+ endif
+ end_unwind_protect
+ else
+ print_usage ();
+ endif
+
+ endfunction
+
+ function [d, f] = myfileparts (x)
+ [d, f, ext] = cellfun (@(x) fileparts (x), x, "UniformOutput", false);
+ f = cellfun (@(x, y) sprintf ("%s%s", x, y), f, ext, ...
+ "UniformOutput", false);
+ idx = cellfun (@(x) isdir (x), x);
+ d(idx) = "";
+ f(idx) = x(idx);
+ endfunction
*** ./scripts/miscellaneous/Makefile.in.orig22 2007-09-05 17:53:19.834645477
+0200
--- ./scripts/miscellaneous/Makefile.in 2007-09-05 22:42:04.087123713 +0200
***************
*** 24,33 ****
compare_versions.m computer.m copyfile.m cputime.m \
delete.m dir.m doc.m dos.m dump_prefs.m \
fileattrib.m fileparts.m flops.m fullfile.m getfield.m gunzip.m \
! inputname.m ispc.m isunix.m license.m list_primes.m ls.m \
ls_command.m menu.m mex.m mexext.m mkoctfile.m movefile.m \
news.m orderfields.m pack.m paren.m parseparams.m \
! run.m semicolon.m setfield.m single.m substruct.m tar.m \
tempdir.m tempname.m texas_lotto.m unix.m unpack.m untar.m \
unzip.m ver.m version.m warning_ids.m xor.m zip.m
--- 24,33 ----
compare_versions.m computer.m copyfile.m cputime.m \
delete.m dir.m doc.m dos.m dump_prefs.m \
fileattrib.m fileparts.m flops.m fullfile.m getfield.m gunzip.m \
! gzip.m inputname.m ispc.m isunix.m license.m list_primes.m ls.m \
ls_command.m menu.m mex.m mexext.m mkoctfile.m movefile.m \
news.m orderfields.m pack.m paren.m parseparams.m \
! run.m semicolon.m setfield.m single.m substruct.m swapbytes.m tar.m \
tempdir.m tempname.m texas_lotto.m unix.m unpack.m untar.m \
unzip.m ver.m version.m warning_ids.m xor.m zip.m
*** ./src/DLD-FUNCTIONS/bsxfun.cc.orig22 2007-09-05 17:45:44.314188192
+0200
--- ./src/DLD-FUNCTIONS/bsxfun.cc 2007-09-05 22:37:55.170799340 +0200
***************
*** 0 ****
--- 1,483 ----
+ /*
+
+ Copyright (C) 2007 David Bateman
+
+ This file is part of Octave.
+
+ Octave is free software; you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation; either version 2, or (at your option) any
+ later version.
+
+ Octave is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Octave; see the file COPYING. If not, write to the Free
+ Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ 02110-1301, USA.
+
+ */
+
+ #ifdef HAVE_CONFIG_H
+ #include <config.h>
+ #endif
+
+ #include <string>
+ #include <vector>
+ #include <list>
+
+ #include "lo-mappers.h"
+
+ #include "oct-map.h"
+ #include "defun-dld.h"
+ #include "parse.h"
+ #include "variables.h"
+ #include "ov-colon.h"
+ #include "unwind-prot.h"
+
+ static bool
+ maybe_update_column (octave_value& Ac, const octave_value& A,
+ const dim_vector& dva, const dim_vector& dvc,
+ octave_idx_type i, octave_value_list &idx)
+ {
+ octave_idx_type nd = dva.length ();
+
+ if (i == 0)
+ {
+ idx(0) = octave_value (':');
+ for (octave_idx_type j = 1; j < nd; j++)
+ {
+ if (dva (j) == 1)
+ idx (j) = octave_value (1);
+ else
+ idx (j) = octave_value ((i % dvc(j)) + 1);
+
+ i = i / dvc (j);
+ }
+
+ Ac = A;
+ Ac = Ac.single_subsref ("(", idx);
+ return true;
+ }
+ else
+ {
+ bool is_changed = false;
+ octave_idx_type k = i;
+ octave_idx_type k1 = i - 1;
+ for (octave_idx_type j = 1; j < nd; j++)
+ {
+ if (dva(j) != 1 && k % dvc (j) != k1 % dvc (j))
+ {
+ idx (j) = octave_value ((k % dvc(j)) + 1);
+ is_changed = true;
+ }
+
+ k = k / dvc (j);
+ k1 = k1 / dvc (j);
+ }
+
+ if (is_changed)
+ {
+ Ac = A;
+ Ac = Ac.single_subsref ("(", idx);
+ return true;
+ }
+ else
+ return false;
+ }
+ }
+
+ static void
+ update_index (octave_value_list& idx, const dim_vector& dv, octave_idx_type i)
+ {
+ octave_idx_type nd = dv.length ();
+
+ if (i == 0)
+ {
+ for (octave_idx_type j = nd - 1; j > 0; j--)
+ idx(j) = octave_value (static_cast<double>(1));
+ idx(0) = octave_value (':');
+ }
+ else
+ {
+ for (octave_idx_type j = 1; j < nd; j++)
+ {
+ idx (j) = octave_value (i % dv (j) + 1);
+ i = i / dv (j);
+ }
+ }
+ }
+
+ static void
+ update_index (Array<int>& idx, const dim_vector& dv, octave_idx_type i)
+ {
+ octave_idx_type nd = dv.length ();
+
+ idx(0) = 0;
+ for (octave_idx_type j = 1; j < nd; j++)
+ {
+ idx (j) = i % dv (j);
+ i = i / dv (j);
+ }
+ }
+
+ DEFUN_DLD (bsxfun, args, nargout,
+ " -*- texinfo -*-\n\
+ @deftypefn {Lodable Function} {} bsxfun (@var{f}, @var{a}, @var{b})\n\
+ Applies a binary function @var{f} element-wise to two matrix arguments\n\
+ @var{a} and @var{b}. The function @var{f} must be capable of accepting\n\
+ two column vector arguments of equal length, or one column vector\n\
+ argument and a scalar.\n\
+ \n\
+ The dimensions of @var{a} and @var{b} must be equal or singleton. The\n\
+ singleton dimensions a the matirces will be expanded to the same\n\
+ dimensioanlity as the other matrix.\n\
+ \n\
+ @seealso{arrayfun, cellfun}\n\
+ @end deftypefn")
+ {
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 3)
+ print_usage ();
+ else
+ {
+ octave_function *func = 0;
+ std::string name;
+ std::string fcn_name;
+
+ if (args(0).is_function_handle () || args(0).is_inline_function ())
+ func = args(0).function_value ();
+ else if (args(0).is_string ())
+ {
+ name = args(0).string_value ();
+ fcn_name = unique_symbol_name ("__bsxfun_fcn_");
+ std::string fname = "function y = ";
+ fname.append (fcn_name);
+ fname.append ("(x) y = ");
+ func = extract_function (args(0), "bsxfun", fcn_name, fname,
+ "; endfunction");
+ }
+ else
+ error ("bsxfun: first argument must be a string or function handle");
+
+ if (! error_state)
+ {
+ const octave_value A = args (1);
+ dim_vector dva = A.dims ();
+ octave_idx_type nda = dva.length ();
+ const octave_value B = args (2);
+ dim_vector dvb = B.dims ();
+ octave_idx_type ndb = dvb.length ();
+ octave_idx_type nd = nda;
+
+ if (nda > ndb)
+ dvb.resize (nda, 1);
+ else if (nda < ndb)
+ {
+ dva.resize (ndb, 1);
+ nd = ndb;
+ }
+
+ for (octave_idx_type i = 0; i < nd; i++)
+ if (dva (i) != dvb (i) && dva (i) != 1 && dvb (i) != 1)
+ {
+ error ("bsxfun: dimensions don't match");
+ break;
+ }
+
+ if (!error_state)
+ {
+ // Find the size of the output
+ dim_vector dvc;
+ dvc.resize (nd);
+
+ for (octave_idx_type i = 0; i < nd; i++)
+ dvc (i) = (dva (i) < 1 ? dva (i) : (dvb (i) < 1 ? dvb (i) :
+ (dva (i) > dvb (i) ? dva (i) : dvb (i))));
+
+ if (dva == dvb || dva.numel () == 1 || dvb.numel () == 1)
+ {
+ octave_value_list inputs;
+ inputs (0) = A;
+ inputs (1) = B;
+ retval = feval (func, inputs, 1);
+ }
+ else if (dvc.numel () < 1)
+ {
+ octave_value_list inputs;
+ inputs (0) = A.resize (dvc);
+ inputs (1) = B.resize (dvc);
+ retval = feval (func, inputs, 1);
+ }
+ else
+ {
+ octave_idx_type ncount = 1;
+ for (octave_idx_type i = 1; i < nd; i++)
+ ncount *= dvc (i);
+
+ #define BSXDEF(T) \
+ T result_ ## T; \
+ bool have_ ## T = false;
+
+ BSXDEF(NDArray);
+ BSXDEF(ComplexNDArray);
+ BSXDEF(boolNDArray);
+ BSXDEF(int8NDArray);
+ BSXDEF(int16NDArray);
+ BSXDEF(int32NDArray);
+ BSXDEF(int64NDArray);
+ BSXDEF(uint8NDArray);
+ BSXDEF(uint16NDArray);
+ BSXDEF(uint32NDArray);
+ BSXDEF(uint64NDArray);
+
+ octave_value Ac ;
+ octave_value_list idxA;
+ octave_value Bc;
+ octave_value_list idxB;
+ octave_value C;
+ octave_value_list inputs;
+ Array<int> ra_idx (dvc.length(), 0);
+
+
+ for (octave_idx_type i = 0; i < ncount; i++)
+ {
+ if (maybe_update_column (Ac, A, dva, dvc, i, idxA))
+ inputs (0) = Ac;
+
+ if (maybe_update_column (Bc, B, dvb, dvc, i, idxB))
+ inputs (1) = Bc;
+
+ octave_value_list tmp = feval (func, inputs, 1);
+
+ if (error_state)
+ break;
+
+ #define BSXINIT(T, CLS, EXTRACTOR) \
+ (result_type == CLS) \
+ { \
+ have_ ## T = true; \
+ result_ ## T = \
+ tmp (0). EXTRACTOR ## _array_value (); \
+ result_ ## T .resize (dvc); \
+ }
+
+ if (i == 0)
+ {
+ if (! tmp(0).is_sparse_type ())
+ {
+ std::string result_type = tmp(0).class_name ();
+ if (result_type == "double")
+ {
+ if (tmp(0).is_real_type ())
+ {
+ have_NDArray = true;
+ result_NDArray = tmp(0).array_value ();
+ result_NDArray.resize (dvc);
+ }
+ else
+ {
+ have_ComplexNDArray = true;
+ result_ComplexNDArray =
+ tmp(0).complex_array_value ();
+ result_ComplexNDArray.resize (dvc);
+ }
+ }
+ else if BSXINIT(boolNDArray, "logical", bool)
+ else if BSXINIT(int8NDArray, "int8", int8)
+ else if BSXINIT(int16NDArray, "int16", int16)
+ else if BSXINIT(int32NDArray, "int32", int32)
+ else if BSXINIT(int64NDArray, "int64", int64)
+ else if BSXINIT(uint8NDArray, "uint8", uint8)
+ else if BSXINIT(uint16NDArray, "uint16", uint16)
+ else if BSXINIT(uint32NDArray, "uint32", uint32)
+ else if BSXINIT(uint64NDArray, "uint64", uint64)
+ else
+ {
+ C = tmp (0);
+ C = C.resize (dvc);
+ }
+ }
+ }
+ else
+ {
+ update_index (ra_idx, dvc, i);
+
+ if (have_NDArray)
+ {
+ if (tmp(0).class_name () != "double")
+ {
+ have_NDArray = false;
+ C = result_NDArray;
+ C = do_cat_op (C, tmp(0), ra_idx);
+ }
+ else if (tmp(0).is_real_type ())
+ result_NDArray.insert (tmp(0).array_value(),
+ ra_idx);
+ else
+ {
+ result_ComplexNDArray =
+ ComplexNDArray (result_NDArray);
+ result_ComplexNDArray.insert
+ (tmp(0).complex_array_value(), ra_idx);
+ have_NDArray = false;
+ have_ComplexNDArray = true;
+ }
+ }
+
+ #define BSXLOOP(T, CLS, EXTRACTOR) \
+ (have_ ## T) \
+ { \
+ if (tmp (0).class_name () != CLS) \
+ { \
+ have_ ## T = false; \
+ C = result_ ## T; \
+ C = do_cat_op (C, tmp (0), ra_idx); \
+ } \
+ else \
+ result_ ## T .insert \
+ (tmp(0). EXTRACTOR ## _array_value (), \
+ ra_idx); \
+ }
+
+ else if BSXLOOP(ComplexNDArray, "double", complex)
+ else if BSXLOOP(boolNDArray, "logical", bool)
+ else if BSXLOOP(int8NDArray, "int8", int8)
+ else if BSXLOOP(int16NDArray, "int16", int16)
+ else if BSXLOOP(int32NDArray, "int32", int32)
+ else if BSXLOOP(int64NDArray, "int64", int64)
+ else if BSXLOOP(uint8NDArray, "uint8", uint8)
+ else if BSXLOOP(uint16NDArray, "uint16", uint16)
+ else if BSXLOOP(uint32NDArray, "uint32", uint32)
+ else if BSXLOOP(uint64NDArray, "uint64", uint64)
+ else
+ C = do_cat_op (C, tmp(0), ra_idx);
+ }
+ }
+
+ #define BSXEND(T) \
+ (have_ ## T) \
+ retval (0) = result_ ## T;
+
+ if BSXEND(NDArray)
+ else if BSXEND(ComplexNDArray)
+ else if BSXEND(boolNDArray)
+ else if BSXEND(int8NDArray)
+ else if BSXEND(int16NDArray)
+ else if BSXEND(int32NDArray)
+ else if BSXEND(int64NDArray)
+ else if BSXEND(uint8NDArray)
+ else if BSXEND(uint16NDArray)
+ else if BSXEND(uint32NDArray)
+ else if BSXEND(uint64NDArray)
+ else
+ retval(0) = C;
+ }
+ }
+ }
+
+ if (! fcn_name.empty ())
+ clear_function (fcn_name);
+ }
+
+ return retval;
+ }
+
+ /*
+
+ %!shared a, b, c, f
+ %! a = randn (4, 4);
+ %! b = mean (a, 1);
+ %! c = mean (a, 2);
+ %! f = @minus;
+ %!error(bsxfun (f));
+ %!error(bsxfun (f, a));
+ %!error(bsxfun (a, b));
+ %!error(bsxfun (a, b, c));
+ %!error(bsxfun (f, a, b, c));
+ %!error(bsxfun (f, ones(4, 0), ones(4, 4)))
+ %!assert(bsxfun (f, ones(4, 0), ones(4, 1)), zeros(4, 0));
+ %!assert(bsxfun (f, ones(1, 4), ones(4, 1)), zeros(4, 4));
+ %!assert(bsxfun (f, a, b), a - repmat(b, 4, 1));
+ %!assert(bsxfun (f, a, c), a - repmat(c, 1, 4));
+ %!assert(bsxfun ("minus", ones(1, 4), ones(4, 1)), zeros(4, 4));
+
+ %!shared a, b, c, f
+ %! a = randn (4, 4);
+ %! a(1) *= 1i;
+ %! b = mean (a, 1);
+ %! c = mean (a, 2);
+ %! f = @minus;
+ %!error(bsxfun (f));
+ %!error(bsxfun (f, a));
+ %!error(bsxfun (a, b));
+ %!error(bsxfun (a, b, c));
+ %!error(bsxfun (f, a, b, c));
+ %!error(bsxfun (f, ones(4, 0), ones(4, 4)))
+ %!assert(bsxfun (f, ones(4, 0), ones(4, 1)), zeros(4, 0));
+ %!assert(bsxfun (f, ones(1, 4), ones(4, 1)), zeros(4, 4));
+ %!assert(bsxfun (f, a, b), a - repmat(b, 4, 1));
+ %!assert(bsxfun (f, a, c), a - repmat(c, 1, 4));
+ %!assert(bsxfun ("minus", ones(1, 4), ones(4, 1)), zeros(4, 4));
+
+ %!shared a, b, c, f
+ %! a = randn (4, 4);
+ %! a(end) *= 1i;
+ %! b = mean (a, 1);
+ %! c = mean (a, 2);
+ %! f = @minus;
+ %!error(bsxfun (f));
+ %!error(bsxfun (f, a));
+ %!error(bsxfun (a, b));
+ %!error(bsxfun (a, b, c));
+ %!error(bsxfun (f, a, b, c));
+ %!error(bsxfun (f, ones(4, 0), ones(4, 4)))
+ %!assert(bsxfun (f, ones(4, 0), ones(4, 1)), zeros(4, 0));
+ %!assert(bsxfun (f, ones(1, 4), ones(4, 1)), zeros(4, 4));
+ %!assert(bsxfun (f, a, b), a - repmat(b, 4, 1));
+ %!assert(bsxfun (f, a, c), a - repmat(c, 1, 4));
+ %!assert(bsxfun ("minus", ones(1, 4), ones(4, 1)), zeros(4, 4));
+
+ %!shared a, b, c, f
+ %! a = randn (4, 4);
+ %! b = a (1, :);
+ %! c = a (:, 1);
+ %! f = @(x, y) x == y;
+ %!error(bsxfun (f));
+ %!error(bsxfun (f, a));
+ %!error(bsxfun (a, b));
+ %!error(bsxfun (a, b, c));
+ %!error(bsxfun (f, a, b, c));
+ %!error(bsxfun (f, ones(4, 0), ones(4, 4)))
+ %!assert(bsxfun (f, ones(4, 0), ones(4, 1)), zeros(4, 0, "logical"));
+ %!assert(bsxfun (f, ones(1, 4), ones(4, 1)), ones(4, 4, "logical"));
+ %!assert(bsxfun (f, a, b), a == repmat(b, 4, 1));
+ %!assert(bsxfun (f, a, c), a == repmat(c, 1, 4));
+
+ %!shared a, b, c, d, f
+ %! a = randn (4, 4, 4);
+ %! b = mean (a, 1);
+ %! c = mean (a, 2);
+ %! d = mean (a, 3);
+ %! f = @minus;
+ %!error(bsxfun (f, ones([4, 0, 4]), ones([4, 4, 4])));
+ %!assert(bsxfun (f, ones([4, 0, 4]), ones([4, 1, 4])), zeros([4, 0, 4]));
+ %!assert(bsxfun (f, ones([4, 4, 0]), ones([4, 1, 1])), zeros([4, 4, 0]));
+ %!assert(bsxfun (f, ones([1, 4, 4]), ones([4, 1, 4])), zeros([4, 4, 4]));
+ %!assert(bsxfun (f, ones([4, 4, 1]), ones([4, 1, 4])), zeros([4, 4, 4]));
+ %!assert(bsxfun (f, ones([4, 1, 4]), ones([1, 4, 4])), zeros([4, 4, 4]));
+ %!assert(bsxfun (f, ones([4, 1, 4]), ones([1, 4, 1])), zeros([4, 4, 4]));
+ %!assert(bsxfun (f, a, b), a - repmat(b, [4, 1, 1]));
+ %!assert(bsxfun (f, a, c), a - repmat(c, [1, 4, 1]));
+ %!assert(bsxfun (f, a, d), a - repmat(d, [1, 1, 4]));
+ %!assert(bsxfun ("minus", ones([4, 0, 4]), ones([4, 1, 4])), zeros([4, 0,
4]));
+
+ %% The below is a very hard case to treat
+ %!assert(bsxfun (f, ones([4, 1, 4, 1]), ones([1, 4, 1, 4])), zeros([4, 4, 4,
4]));
+
+ */
*** ./src/DLD-FUNCTIONS/typecast.cc.orig22 2007-09-05 17:45:50.767869321
+0200
--- ./src/DLD-FUNCTIONS/typecast.cc 2007-09-05 18:05:43.729801202 +0200
***************
*** 0 ****
--- 1,211 ----
+ /*
+
+ Copyright (C) 2007 David Bateman
+
+ This file is part of Octave.
+
+ Octave is free software; you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation; either version 2, or (at your option) any
+ later version.
+
+ Octave is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with Octave; see the file COPYING. If not, write to the Free
+ Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
+ 02110-1301, USA.
+
+ */
+
+ #ifdef HAVE_CONFIG_H
+ #include <config.h>
+ #endif
+
+ #include <string>
+
+ #include "oct.h"
+
+ template <class LT, class RT>
+ static void
+ typecast (const Array<RT>& x, Array<LT>& y)
+ {
+ octave_idx_type n = x.length ();
+ size_t ns = sizeof (RT);
+ size_t ms = sizeof (LT);
+
+ if (n * ns % ms != 0)
+ error ("typecast: incorrect number of input values to make output
value");
+ else
+ {
+ octave_idx_type m = n * ns / ms;
+ dim_vector dv = x.dims ();
+ for (octave_idx_type i = 0; i < dv.length(); i++)
+ if (dv(i) == n)
+ {
+ dv(i) = m;
+ break;
+ }
+ y.resize (dv);
+ const unsigned char *xp = reinterpret_cast<const unsigned char *>
+ (x.fortran_vec ());
+ unsigned char *yp = reinterpret_cast<unsigned char *>(y.fortran_vec ());
+ for (octave_idx_type i = 0;
+ i < n * static_cast<octave_idx_type>(ns); i++)
+ *yp++ = *xp++;
+ }
+ }
+
+ template <class T>
+ static octave_value
+ typecast (const T& x, std::string type)
+ {
+ octave_value retval;
+ if (type == "uint8")
+ {
+ uint8NDArray y;
+ typecast (x, y);
+ retval = octave_value (y);
+ }
+ else if (type == "uint16")
+ {
+ uint16NDArray y;
+ typecast (x, y);
+ retval = octave_value (y);
+ }
+ else if (type == "uint32")
+ {
+ uint32NDArray y;
+ typecast (x, y);
+ retval = octave_value (y);
+ }
+ else if (type == "uint64")
+ {
+ uint64NDArray y;
+ typecast (x, y);
+ retval = octave_value (y);
+ }
+ else if (type == "int8")
+ {
+ int8NDArray y;
+ typecast (x, y);
+ retval = octave_value (y);
+ }
+ else if (type == "int16")
+ {
+ int16NDArray y;
+ typecast (x, y);
+ retval = octave_value (y);
+ }
+ else if (type == "int32")
+ {
+ int32NDArray y;
+ typecast (x, y);
+ retval = octave_value (y);
+ }
+ else if (type == "int64")
+ {
+ int64NDArray y;
+ typecast (x, y);
+ retval = octave_value (y);
+ }
+ else
+ {
+ NDArray y;
+ typecast (x, y);
+ retval = octave_value (y);
+ }
+
+ return retval;
+ }
+
+ DEFUN_DLD (typecast, args, ,
+ " -*- texinfo -*-\n\
+ @deftypefn {Loadable Function} {} typecast (@var{x}, @var{type})\n\
+ Converts from one datatype to another without changing the underlying\n\
+ data. The argument @var{type} defines the type of the return argument\n\
+ and must be one of 'uint8', 'uint16', 'uint32', 'uint64', 'int8', 'int16',\n\
+ 'int32', 'int64', 'single' or 'double'.\n\
+ \n\
+ An example of the use of typecast on a little-endian machine is\n\
+ \n\
+ @example\n\
+ @group\n\
+ @var{x} = uint16 ([1, 65535]);\n\
+ typecast (@var{x}, 'uint8')\n\
+ @result{} [ 0, 1, 255, 255]\n\
+ @end group\n\
+ @end example\n\
+ @seealso{cast, swapbytes}\n\
+ @end deftypefn")
+ {
+ int nargin = args.length ();
+ octave_value retval;
+
+ if (nargin != 2)
+ print_usage ();
+ else
+ {
+ std::string type = args (1).string_value ();
+
+ if (! error_state)
+ {
+ std::transform (type.begin (), type.end (), type.begin (), tolower);
+
+ if (type == "single")
+ error ("typecast: type 'single' is not supported");
+ else if (type != "uint8" && type != "uint16" &&
+ type != "uint32" && type != "uint64" &&
+ type != "int8" && type != "int16" &&
+ type != "int32" && type != "int64" &&
+ type != "single" && type != "double")
+ error ("typecast: unrecognized or invalid type");
+ else if (args(0).is_sparse_type () || args(0).is_complex_type ())
+ error ("typecast: sparse and complex types are invalid");
+ else
+ {
+ dim_vector dv = args(0).dims ();
+ bool seen = false;
+
+ for (octave_idx_type i = 0; i < dv.length(); i++)
+ if (dv (i) != 1)
+ {
+ if (seen)
+ {
+ error ("typecast: argument must be a vector");
+ break;
+ }
+ else
+ seen = true;
+ }
+
+ if (!error_state)
+ {
+ if (args(0).is_uint8_type ())
+ retval = typecast (args(0).uint8_array_value (), type);
+ else if (args(0).is_uint16_type ())
+ retval = typecast (args(0).uint16_array_value (), type);
+ else if (args(0).is_uint32_type ())
+ retval = typecast (args(0).uint32_array_value (), type);
+ else if (args(0).is_uint64_type ())
+ retval = typecast (args(0).uint64_array_value (), type);
+ else if (args(0).is_int8_type ())
+ retval = typecast (args(0).int8_array_value (), type);
+ else if (args(0).is_int16_type ())
+ retval = typecast (args(0).int16_array_value (), type);
+ else if (args(0).is_int32_type ())
+ retval = typecast (args(0).int32_array_value (), type);
+ else if (args(0).is_int64_type ())
+ retval = typecast (args(0).int64_array_value (), type);
+ else
+ retval = typecast (args(0).array_value (), type);
+ }
+ }
+ }
+ }
+
+ return retval;
+ }
*** ./src/Makefile.in.orig22 2007-09-05 17:46:09.057965542 +0200
--- ./src/Makefile.in 2007-09-05 17:47:33.850774054 +0200
***************
*** 47,53 ****
OPT_HANDLERS := DASPK-opts.cc DASRT-opts.cc DASSL-opts.cc \
LSODE-opts.cc NLEqn-opts.cc Quad-opts.cc
! DLD_XSRC := balance.cc besselj.cc betainc.cc cellfun.cc chol.cc \
ccolamd.cc colamd.cc colloc.cc conv2.cc convhulln.cc daspk.cc \
dasrt.cc dassl.cc det.cc dispatch.cc eig.cc expm.cc \
fft.cc fft2.cc fftn.cc fftw.cc filter.cc find.cc fsolve.cc \
--- 47,53 ----
OPT_HANDLERS := DASPK-opts.cc DASRT-opts.cc DASSL-opts.cc \
LSODE-opts.cc NLEqn-opts.cc Quad-opts.cc
! DLD_XSRC := balance.cc besselj.cc betainc.cc bsxfun.cc cellfun.cc chol.cc \
ccolamd.cc colamd.cc colloc.cc conv2.cc convhulln.cc daspk.cc \
dasrt.cc dassl.cc det.cc dispatch.cc eig.cc expm.cc \
fft.cc fft2.cc fftn.cc fftw.cc filter.cc find.cc fsolve.cc \
***************
*** 56,65 ****
lu.cc luinc.cc matrix_type.cc md5sum.cc minmax.cc pinv.cc qr.cc \
quad.cc qz.cc rand.cc regexp.cc schur.cc sort.cc sparse.cc \
spchol.cc spdet.cc spfind.cc spkron.cc splu.cc spparms.cc spqr.cc \
! sqrtm.cc svd.cc syl.cc symrcm.cc time.cc tsearch.cc urlwrite.cc \
! __contourc__.cc __delaunayn__.cc __dsearchn__.cc __gnuplot_raw__.l \
! __glpk__.cc __lin_interpn__.cc __pchip_deriv__.cc __qp__.cc \
! __voronoi__.cc
DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))
--- 56,65 ----
lu.cc luinc.cc matrix_type.cc md5sum.cc minmax.cc pinv.cc qr.cc \
quad.cc qz.cc rand.cc regexp.cc schur.cc sort.cc sparse.cc \
spchol.cc spdet.cc spfind.cc spkron.cc splu.cc spparms.cc spqr.cc \
! sqrtm.cc svd.cc syl.cc symrcm.cc time.cc tsearch.cc typecast.cc \
! urlwrite.cc __contourc__.cc __delaunayn__.cc __dsearchn__.cc \
! __gnuplot_raw__.l __glpk__.cc __lin_interpn__.cc __pchip_deriv__.cc \
! __qp__.cc __voronoi__.cc
DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))
2007-09-06 David Bateman <address@hidden>
* interpreter/system.m: Document gzip.
* interpreter/container.txi: Document celldisp.
* interpreter/matrix.txi: Document bsxfun.
* interpreter/data.txi: Document typecast and swapbytes.
2007-09-06 David Bateman <address@hidden>
* Array-util.cc (increment_index): dimensions can have singleton
trailing dimensions.
2007-09-06 David Bateman <address@hidden>
* general/celldisp.m: New function.
* general/Makefile.in (SOURCES): Add celldisp.m.
* miscellaneous/swapbytes.m: New function.
* miscellaneous/gzip.m: New function.
* miscellaneous/Makefile.in (SOURCES): Add swapbytes.m and gzip.m.
2007-09-06 David Bateman <address@hidden>
* DLD-FUNCTIONS/bsxfun.cc: New function.
* DLD-FUNCTIONS/typecast.cc: New function.
* Makefile.in (DLD_XSRC): Add bsxfun.cc and typecast.cc.
- Five other functions that are in Matlab core ported to Octave,
David Bateman <=