Octave, GNU Octave, Matlab, Scientific Computing, Language, Interpreter, Compiler, C++, LAPACK, Fortran, Fun , GNU

Showing posts with label Paul Kienzle. Show all posts
Showing posts with label Paul Kienzle. Show all posts

Friday, January 18, 2008

Octave-2-Matlab and back (2)

Like I mentioned earlier here. I finally managed to put together a AST-visitor pattern that
loads the user-script parses it, and gives it to the conversion tree-walker code. This implements
the changes mentioned in the oct2mat package (find it here) by Paul Kienzle .

The email announcing this code and patch was posted

Hello there,

I have a oct2mat like the one Paul Kienzle mentions in his
README file based on the AST walker. Utilizing Octave's
parse tree to do the conversion is very helpful and accurate
to a larger degree. Its still a work in progress (my disclaimer).

Large parts of the code are directly derived from Octave
source written by JWE.

This script is essentially to help deploy code to your colleagues
who want to use it on Matlab systems for whatever reason. Primary
development is still expected to be based on an Octave based
environment.

You can invoke the script as,

octave:1> oct2mat("full/path/my_cool_octave_function.m",1)
octave:2> oct2mat("full/path/my_cool_octave_script.m")

and the code is output to stdout, where from you can use some shell
scripts to redirect the output to another matlab directory.

I have attached the files. Im waiting for JWE to accept some
global export patches to access some parser information. Once that is
done, I plan to add it to the extras branch of octave-forge.

This code works only on the cvs version of Octave, as a external add-on.

Limitations include

1. Can treat only function files or script files not both.
2. Pollutes your command history, as it pipes out to octave_stdout.
3. Cannot handle a script file with function definitions
4. Doesnot write to a separate directory or file, just all output on
stdout. So you need to create a driver program.
5. Not sure if it works correctly for all cases.

Best,
-Muthu



The Code / Tree Walker

I also made patches against existing code from Octave tree (maybe obsolete right now), so
as to help integrate my tree walker into the existing visitor.

Index: pt-pr-code.cc
===================================================================
RCS file: /cvs/octave/src/pt-pr-code.cc,v
retrieving revision 1.54
diff -u -p -r1.54 pt-pr-code.cc
--- pt-pr-code.cc 12 Oct 2007 21:27:33 -0000 1.54
+++ pt-pr-code.cc 10 Jan 2008 05:56:20 -0000
@@ -35,6 +35,240 @@ along with Octave; see the file COPYING.
#include "pr-output.h"
#include "pt-all.h"

+
+static std::string
+matlab_unary_op_as_string (octave_value::unary_op op)
+{
+ std::string retval;
+
+ switch (op)
+ {
+ case octave_value::op_not:
+ retval = "~";
+ break;
+
+ case octave_value::op_uplus:
+ retval = "+";
+ break;
+
+ case octave_value::op_uminus:
+ retval = "-";
+ break;
+
+ case octave_value::op_transpose:
+ retval = ".'";
+ break;
+
+ case octave_value::op_hermitian:
+ retval = "'";
+ break;
+
+ case octave_value::op_incr:
+ retval = "++";
+ break;
+
+ case octave_value::op_decr:
+ retval = "--";
+ break;
+
+ default:
+ retval = "
";
+ }
+
+ return retval;
+}
+
+static std::string
+matlab_binary_op_as_string (octave_value::binary_op op)
+{
+ std::string retval;
+
+ switch (op)
+ {
+ case octave_value::op_add:
+ retval = "+";
+ break;
+
+ case octave_value::op_sub:
+ retval = "-";
+ break;
+
+ case octave_value::op_mul:
+ retval = "*";
+ break;
+
+ case octave_value::op_div:
+ retval = "/";
+ break;
+
+ case octave_value::op_pow:
+ retval = "^";
+ break;
+
+ case octave_value::op_ldiv:
+ retval = "\\";
+ break;
+
+ case octave_value::op_lshift:
+ retval = "<<"; + break; + + case octave_value::op_rshift: + retval = ">>";
+ break;
+
+ case octave_value::op_lt:
+ retval = "<"; + break; + + case octave_value::op_le: + retval = "<="; + break; + + case octave_value::op_eq: + retval = "=="; + break; + + case octave_value::op_ge: + retval = ">=";
+ break;
+
+ case octave_value::op_gt:
+ retval = ">";
+ break;
+
+ case octave_value::op_ne:
+ retval = "~=";
+ break;
+
+ case octave_value::op_el_mul:
+ retval = ".*";
+ break;
+
+ case octave_value::op_el_div:
+ retval = "./";
+ break;
+
+ case octave_value::op_el_pow:
+ retval = ".^";
+ break;
+
+ case octave_value::op_el_ldiv:
+ retval = ".\\";
+ break;
+
+ case octave_value::op_el_and:
+ retval = "&";
+ break;
+
+ case octave_value::op_el_or:
+ retval = "|";
+ break;
+
+ case octave_value::op_struct_ref:
+ retval = ".";
+ break;
+
+ default:
+ retval = "
";
+ }
+
+ return retval;
+}
+
+static void
+do_oct2mat_oper(std::ostream &os, tree_binary_expression &expr)
+{
+ std::string oprval="
";
+
+ // this converts code into whatever form
+ // oct2mat requires.
+ if ( expr.is_boolean_expression() )
+ {
+ const tree_boolean_expression & tbr
+ = dynamic_cast
(expr);
+
+ oprval="
";
+
+ switch ( tbr.op_type() )
+ {
+ case tree_boolean_expression::bool_and:
+ oprval = "&";
+ break;
+
+ case tree_boolean_expression::bool_or:
+ oprval = "|";
+ break;
+
+ defualt:
+ oprval = "
";
+ }
+ }
+ else
+ {
+ oprval=matlab_binary_op_as_string(expr.op_type());
+ }
+
+ os << " " << oprval="">";
+
+ oprval=matlab_unary_op_as_string(expr.op_type());
+ os << " " <<>accept(omat);
+ buf_str = buf.str();
+ os << idx="0;">= 0);
+
+ s << ""; + if (printing_newlines) + { + if (beginning_of_line) + { + s << i =" 0;" beginning_of_line =" false;">accept (*this);

- os << " " << op2 =" expr.rhs">accept (*this);

@@ -237,7 +474,10 @@ tree_print_code::visit_simple_for_comman

indent ();

- os << "endfor"; + if ( convert_to_matlab ) + os << "end"; + else + os << "endfor"; } void @@ -280,7 +520,10 @@ tree_print_code::visit_complex_for_comma indent (); - os << "endfor"; + if ( convert_to_matlab ) + os << "end"; + else + os << "endfor"; } void @@ -401,7 +644,10 @@ tree_print_code::visit_octave_user_funct indent (); - os << "endfunction"; + if ( convert_to_matlab ) + os << "%% endfunction"; + else + os << "endfunction"; newline (); } @@ -414,7 +660,17 @@ tree_print_code::visit_identifier (tree_ print_parens (id, "("); std::string nm = id.name (); - os << (nm.empty () ? std::string ("(empty)") : nm); + + if ( convert_to_matlab ) + { + std::map
::iterator iter= identifier_lut.find(nm);
+ if ( iter != identifier_lut.end() )
+ os <<>second ;
+ else
+ os << (nm.empty () ? std::string ("(empty)") : nm); + } + else + os << (nm.empty () ? std::string ("(empty)") : nm); print_parens (id, ")"); } @@ -459,7 +715,10 @@ tree_print_code::visit_if_command (tree_ indent (); - os << "endif"; + if ( convert_to_matlab ) + os << "end"; + else + os << "endif"; } void @@ -693,11 +952,75 @@ tree_print_code::visit_no_op_command (tr void tree_print_code::visit_constant (tree_constant& val) { + std::string buf_str; + std::ostringstream buf; + int idx = 0; + octave_value sval=val.rvalue(); indent (); print_parens (val, "("); - val.print_raw (os, true, print_original_text); + + if ( convert_to_matlab ) + { + if ( sval.is_string() ) + { + buf_str=sval.string_value(); + // see if you need to escape the string stxt. + os << "'"; + for( idx=0; idx < buf_str =" buf.str();" idx="0;">accept (*this);

- os << e =" expr.operand" s =" indent_to_string">accept (*this);
+ // print a commented code.
+ if ( convert_to_matlab )
+ {
+ std::string ind_prfx = indent_to_string( );
+ do_print_commented_code(unwind_protect_code,os,ind_prfx);
+ }
+ else
+ unwind_protect_code->accept (*this);

decrement_indent_level ();
}
@@ -1022,7 +1374,10 @@ tree_print_code::visit_unwind_protect_co

indent ();

- os << "unwind_protect_cleanup"; + if ( convert_to_matlab ) + os << "%%unwind_protect_cleanup"; + else + os << "unwind_protect_cleanup"; newline (); @@ -1032,7 +1387,13 @@ tree_print_code::visit_unwind_protect_co { increment_indent_level (); - cleanup_code->accept (*this);
+ if ( convert_to_matlab )
+ {
+ std::string ind_prfx = indent_to_string( );
+ do_print_commented_code(cleanup_code,os,ind_prfx);
+ }
+ else
+ cleanup_code->accept (*this);

decrement_indent_level ();
}
@@ -1041,7 +1402,10 @@ tree_print_code::visit_unwind_protect_co

indent ();

- os << "end_unwind_protect"; + if ( convert_to_matlab ) + os << "%%end_unwind_protect"; + else + os << "end_unwind_protect"; } void @@ -1075,7 +1439,10 @@ tree_print_code::visit_while_command (tr indent (); - os << "endwhile"; + if ( convert_to_matlab ) + os << "end"; + else + os << "endwhile"; } void @@ -1192,7 +1559,10 @@ tree_print_code::print_comment_elt (cons if (c == '\n') { if (prev_char_was_newline) - os << "##"; + if ( convert_to_matlab ) + os << "%" ; + else + os << "##"; newline (); @@ -1206,7 +1576,11 @@ tree_print_code::print_comment_elt (cons indent (); - os << "##"; + + if ( convert_to_matlab ) + os << "%%"; + else + os << "##"; if (! (isspace (c) || c == '!')) os << " "; Index: pt-pr-code.h =================================================================== RCS file: /cvs/octave/src/pt-pr-code.h,v retrieving revision 1.36 diff -u -p -r1.36 pt-pr-code.h --- pt-pr-code.h 12 Oct 2007 21:27:33 -0000 1.36 +++ pt-pr-code.h 10 Jan 2008 05:56:20 -0000 @@ -26,6 +26,7 @@ along with Octave; see the file COPYING. #include

#include

+#include


#include "comment-list.h"
#include "pt-walk.h"
@@ -41,14 +42,23 @@ public:

tree_print_code (std::ostream& os_arg,
const std::string& pfx = std::string (),
- bool pr_orig_txt = true)
+ bool pr_orig_txt = true, bool to_matlab = false)
: os (os_arg), prefix (pfx), nesting (),
- print_original_text (pr_orig_txt),
+ print_original_text (pr_orig_txt), convert_to_matlab( to_matlab ),
curr_print_indent_level (0), beginning_of_line (true),
printing_newlines (true)
{
// For "none".
nesting.push ('n');
+ if ( convert_to_matlab )
+ {
+ identifier_lut["SEEK_CUR"]="0";
+ identifier_lut["SEEK_END"]="1";
+ identifier_lut["SEEK_SET"]="-1";
+ identifier_lut["usage"]="error";
+ identifier_lut["print_usage"]="error";
+ identifier_lut["__error_text__"]="lasterr";
+ }
}

~tree_print_code (void) { }
@@ -147,6 +157,12 @@ private:

bool print_original_text;

+ // flag controls if tree print code converts to Matlab.
+ bool convert_to_matlab;
+
+ //
+ std::map
identifier_lut;
+
// Current indentation.
int curr_print_indent_level;

@@ -176,6 +192,8 @@ private:

void print_indented_comment (octave_comment_list *comment_list);

+ std::string indent_to_string( );
+
// Must create with an output stream!

tree_print_code (void);




I still have bugs in the driver routine, but hey it works!.
-Muthu

Tuesday, June 26, 2007

Speeding up / Vectorization of Octave code

There was a really interesting posts on the Help-Octave lists yesterday (25th) June where
a vectorization paradigm was elucidated. Here it is for the non-Octavers...

Paul Kienzle to Rob Teasdale
date Jun 25, 2007 9:34 PM
subject Re: Request help speeding up script


On Jun 25, 2007, at 5:23 PM, Rob Teasdale wrote:
> Hi all,
>
> I am still learning Octave, and I am still discovering new ways of
> optimising code. I have managed to vectorise a lot of my 'for' loops,
> and I have managed to speed my small script up dramatically, although
> I am having trouble with this last code snippet below. What I am
> trying to achieve is to calculate the area of a quadrilateral by the
> vector cross product of its diagonals. X, Y and Z contain the points
> that describe my mesh. I am sure that there is a way to vectorise at
> least part of this, and I would really appreciate any suggestions,
>
> [R C] = size(X); % all matrices the same size
> row = 0;
> for i = 1:(R-1)
> for j = 1:(C-1)
> P4 = [X(i,j) Y(i,j) Z(i,j)];
> P3 = [X(i,j+1) Y(i,j+1) Z(i,j+1)];
> P2 = [X(i+1,j+1) Y(i+1,j+1) Z(i+1,j+1)];
> P1 = [X(i+1,j) Y(i+1,j) Z(i+1,j)];
> p = P2-P4;
> q = P3-P1;
> row++
> r(row,:) = 0.5 * (cross(p,q)); % calculate the area
> endfor
> endfor



Here is a mechanical translation of your loop, producing 3-column
matrices p and q:

[m,n] = size(X);
I = 1:m-1;
J = 1:n-1;
P4 = [X(I,J)(:), Y(I,J)(:), Z(I,J)(:)];
P3 = [X(I,J+1)(:),Y(I,J+1)(:), Z(I,J+1)(:)];
P2 = [X(I+1,J+1)(:),Y(I+1,J+1)(:),Z(I+1,J+1)(:)];
P1 = [X(I+1,J)(:),Y(I+1,J)(:),Z(I+1,J)(:)];
p = P2-P4;
q = P3-P1;
r = 0.5 * cross(p,q);

Because you traverse by column rather than by row the results will be in a different order than you give. If for some reason you need your particular order, you can transpose X,Y,Z and swap P3 and P1.

- Paul


Now I wish I could do that, and what we call optimization.
-Muthu

Creative Commons License