outer product implementation

Muthiah Annamalai <[email protected]> Sun, 10 Sep 2006 12:59:51 -0500
Newsgroups gmane.comp.gnu.octave.sources
Message-ID <1157911192.5512.20.camel@localhost>
--=-FeeHII4f1E7Y0AidcIL8
Content-Type: text/plain
Content-Transfer-Encoding: 7bit

Paul,
I have fixed the 2nd version of the outer product code as you had asked
me
to using C++ style inline functions etc & TexInfo documentation.

The older code can be found at.
http://www.cae.wisc.edu/pipermail/octave-sources/2006-April/000000.html

Im attaching the new code & the same test case.
Thanks
Muthu


--=-FeeHII4f1E7Y0AidcIL8
Content-Disposition: attachment; filename=outer.cpp
Content-Type: text/x-c++src; name=outer.cpp; charset=UTF-8
Content-Transfer-Encoding: 7bit

/* 
 * (C) 2006, September, Muthiah Annamalai. <[email protected]>
 *    An implementation of the 'outer'-product function as specified in the
 *    octave-projects page, at www.octave.org.
 *
 *    This program 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 of the License, or
 *    (at your option) any later version.
 *
 *    This program 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 this program; if not, write to the Free Software
 *    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */


#include<iostream>
#include<octave/oct.h>
#include<octave/parse.h>
#include<octave/dynamic-ld.h>
#include<octave/oct-map.h>
#include<octave/oct-stream.h>
#include<octave/ov.h>
#include<octave/parse.h>

#include<octave/Matrix.h>
#include<octave/ov-cx-mat.h>
#include<octave/ov-list.h>
#include<octave/ov-base-mat.h>

#include<octave/quit.h>
#include<octave/error.h>

inline static bool is_vector(octave_value arg) {
	return ((arg.is_matrix_type()&& !arg.is_char_matrix() &&  (arg.columns()==1 || arg.rows()==1))? true : false);
}


DEFUN_DLD(outer,args,,
"-*- texinfo -*-\n"
"@deftypefn {Loadable Function} {@var{outer_product} =} outer (@var{x},@var{y})\n"
"@deftypefnx {Loadable Function} {@var{outer_product} =} outer (@var{x},@var{y},@var{f})\n"
"\n"
"function outer(VecX,VecY,opt:function)\n"
"Return a @var{M}x@var{N} matrix with @var{f}(@var{x},@var{y}) evaluated at x,y position\n"
"for each @var{x},@var{y} in Vector @var{X},Vector @var{Y} respectively\n"
"If function is not provided, use '*' as default operation\n"
"Function can be provided using argument @var{f} \n"
"using reference operator (@@) only.\n"
"\n"
"example: outer(@var{[1.0 2.0 3.0]},@var{[2.0 3.0 4.0]},@@gcd)\n"
"compute the 3x3 matrix with each element in it being the\n"
"GCD of the corresponding x and y elements of the matrix\n"
"\n"
"@end deftypefn\n"
"@seealso{}")
{
  int m=0,n=0,aLength=0;
  ColumnVector cv[2];
  bool function_present=false;
  octave_function *func=NULL;
  octave_value_list erval;
  const char *func_name=NULL;
  

  bool is_commutative=false;

  aLength=args.length();

  if(aLength < 2)
    {
      error("Usage: outer(VecX,VecY,opt:function)");
      return erval;
    }

  if(!(is_vector(args(0)) || is_vector(args(1))))
    {
      error("Usage: outer(VecX,VecY,opt:function);\n Invalid argument arg- [1 or 2]Expected Row/Column Vector");
      return erval;
    }

  /* 
   * FIXME: [Support for String-Name functions. Only references supported now]
   * Also allow functions to be passed by their names, and then replace the
   * name with a reference to the function-, later.
   * hope to use commutativity of the operators someday
   */
  if(aLength>=3 && (args(2).is_function_handle() /* ||args(2).is_string() */))
    {
      function_present=true;
      is_commutative=false;
      
      func=args(2).function_value();
    }

  cv[0]=ColumnVector(args(0).vector_value());
  cv[1]=ColumnVector(args(1).vector_value());

  m=cv[0].length();
  n=cv[1].length();

  Matrix rmat(m,n);
  octave_value_list ovl;
  
  /*
   * If operator is commutative, then
   * do we want to handle this as a special
   * case?
   *
   */

  if(function_present && func!=NULL)
    {
      /* try a dummy invoke on the function,
       and verify if the return value is there. Otherwise
      we report and error & die.*/

      ovl(0)=0;
      ovl(1)=0;

      func_name = func->name().c_str();

      if(feval(func,ovl).length() < 1)
	{
	  error("The function %s must return atleast one value, on invocation\n",func_name);
	  return octave_value(-1);
	}

      if(error_state)
	{
	  error("Exception: Cannot find function specified");
	  return erval;
	}
      
      /* FIXME: <Vectorize Code>
       * Paul says do this thing as a single loop.
       * Iterate over Columns, taking one full column at a time, 
       * and save the result on the columns of the matrix.
       * Atleast this is possible for builtins that take arguments
       * over columns/rows.
       */
      for(int itr=0;itr<m;itr++)
	for(int icol=0;icol<n;icol++)
	  {

	    OCTAVE_QUIT;
	      
	    ovl(0)=cv[0].elem(itr);
	    ovl(1)=cv[1].elem(icol);

	    rmat.elem(itr,icol)=feval(func,ovl)(0).double_value();
	    if(error_state)
	      {
		error("Exception: while evaluating %s with arguments %g %g\n",func_name,cv[0].elem(itr),cv[1].elem(icol));
		return erval;
	      }
	  }
    }
  else 
    {
      /* do a simple product */
      rmat=cv[0]*cv[1].transpose();
    }

  return octave_value(rmat);
}
/* Make command:
g++ -fpic outer.cpp -shared -o outer.oct -Wall -ggdb -Wall -Wunused -Wconversion \
-fno-exceptions -DDEBUG=1  `mkoctfile -p INCFLAGS` 
*/

--=-FeeHII4f1E7Y0AidcIL8
Content-Disposition: attachment; filename=outertest.m
Content-Type: text/x-octave; name=outertest.m; charset=UTF-8
Content-Transfer-Encoding: 7bit

#! /usr/bin/octave -q
%outer
disp('Help on Outer')
help outer
pause
disp('Product- testing ;-) Imean this * product')
outer([1 2],[1 2 3 4])
outer([1],[1 2 3])
outer([1 2 3],[1])
outer([1 2 3],[1 2 3])
%outer([1],[2]) %error this is caught.

disp('Functions being Invoked')

disp('Min')
outer(rand(1,10),rand(1,10),@min)

disp('IsEqual')
outer([1 2],[1 2 3 4],@isequal)

disp('Max')
outer([1],[1 2 3],@max)

disp('GCD ')
outer([1.0 2.0 3.0],[2.0 3.0 4.0],@gcd)
outer([14 20 35 40],[7],@gcd)


disp('Square Matrix')
outer([1:5],[1:5])

disp('One arg scalar other vector')
outer(3,[1:5])
outer([1:5],3)

disp('functions that dont return an argument! ')
try
  outer([14 20 35 40],[7],@disp) %functions that dont return an argument! 
catch
  disp('Exception')
end

disp('non-existent functions throw errors at Octave prompt itself.');
try
outer([14 20 35 40],[7],@bwand) %non-existent functions throw errors at Octave prompt itself.
catch
  disp('Exception')
end

try
disp('Both Scalar')
outer([1],[2]) %error 
catch
  disp('Exception')
end

--=-FeeHII4f1E7Y0AidcIL8
Content-Type: text/plain; charset="us-ascii"
MIME-Version: 1.0
Content-Transfer-Encoding: 7bit
Content-Disposition: inline

_______________________________________________
Octave-sources mailing list
[email protected]
https://www.cae.wisc.edu/mailman/listinfo/octave-sources

--=-FeeHII4f1E7Y0AidcIL8--