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

Tuesday, September 25, 2007

Octave & Set theory

Octave can handle first-order sets, by which I mean, all sets excluding the 'x' (cross, cartesian) operators. This means, we cannot treat power sets properly.

The various set operators and functions that implement them are:

Operation

Function

set definition

create_set

Union

union

Subset

ismember

X’s complement w.r.t Y (Y-X )

complement(X,Y)

X-Y also written (X\Y)

setdiff(X,Y)

(A-B) U (B-A)

setxor(A,B)


Some interesting features to add, would be the capacity to deal with ordered pairs, i.e the same cross product we discussed earlier. We must also be able to define functions on sets, that do operations on such sets; Also we could have subset function return proper & improper subset differences.

I would also suggest generators, for finite sets, and ability to deal with groups. But then we are moving outside the scope of the set theory; this completely limits moves us to having a group theory implementation in Octave, instead of 'solving' this problem. Maybe I dont understand this very well, anyway. Cardinality operator is just the length() function. I also dont write
here the aim of the function 'intersect' that is more like a generic CS problem, as its not a
defined set theoretic operator.

Anyway, trying to force all of mathematics down into a single linear algebra package, seems like masochism.

Polynomials, Sequences...

Polynomial families that make up power series solutions to standard ('special') differential
equations, are useful in a variety of engineering situations (see: Advanced Engineering Mathematics; Kreyzig, Erwin).

Matlab functions for associated Legendre functions; many more functions are present
in Matlab file exchange, but these are either unlicensed or distributed under uncertain
terms; anycase not widely available.

I dont know of better Matlab functions resource for orthogonal polynomial function sets.
I wanted to say, we need the recursive/iterative implementation of these functions, for
GNU Octave.

-Muthu

Friday, September 21, 2007

Information theory toolbox

I've been waiting to write and integrate some of similar code posted,
into the information-theory toolbox today, and this is going to be
possible sometime soon. It extends 3 functions, to make use of sequence
data without having to send in a list of probabilities to compute entropy,
mutual information, and related variables.

An engineering/physics way of thinking about problems includes calculations
and looking at the mental-"picture"; a la the visualization process, is where I
think tools like Octave, Matlab help. They supplement the visualization.

I'm also thinking about making a bunch of functions for doing 1D-matrix optics
using the ABCD matrix formalism. Infact, just have to plugin most of the equations
from the text right into Octave, and its just a simple exercise. Well, then you
could use it to calculate ray-propagation through a optic system, accurate to
the paraxial approximation. Next, we could also study the plane wave gaussian
intensity distributions from the source object plane to the image plane, with
such a system. Ofcourse we have to pick conventions, which makes it hard
for someone (who uses other conventions).

Simply calculating path of each point in the object plane, onto the image place,
we can calculate the intensity of the input beam at the output plane. I think that
is interesting to get a qualitative picture of many Optical instruments, in 1D.

Anyway, I have more plans on the drawing board, than I have time for possibly.
Well, the profiler, oct2mat translator and debugger seem to be shelved indefnitely.
Then, I wanted a small 2D/1D FDTD solver for Maxwell's equations for Octave.
This would be super-cool in my opinion, and sometimes directly useful for my
research work as well!

Best,
-Muthu

Wednesday, June 27, 2007

Octave Graphics, XML and GUI's : Beginning (yet another... )

Im reposting some stuff from the June, 2007 Octave-Maintainers mailing list here,

HOWTO add a new graphics property

Shai Ayal shaiay at gmail.com
Sat Jun 16 22:41:48 CDT 2007

Hi all,

I wrote this up just to have some simple guide with all the necessary
steps on how to add a new graphics property. Please have a look and
correct it if necessary. Hopefully when finalized this howto will
enable more people to contribute code.

Shai

how to add property pppp of type tttt to axes:
all functions mentioned below are in src/graphics.cc and src/graphics.h

in axes_properties class definition, add:
OCTAVE_GRAPHICS_PROPERTY (tttt, pppp);

axes::axes_properties::axes_properties
add code to initialize

axes::axes_properties::set
add code to set ppp from an octave_value

axes::axes_properties::set_defaults
add code to set the default value

axes::axes_properties::get (void) const
add code to return as a cell

axes::axes_properties::get (const property_name& name) const
add code to return as an octave_value

property_list::pval_map_type axes::axes_properties::factory_defaults
add the factory default vlaue for pppp

if the property is a child object, also modify

void axes::axes_properties::delete_children
void axes::axes_properties::remove_child

Also, I noticed this last month and a littler earlier to that, many new projects on Octave, including
  1. Qt-Octave, developed by from Spain
  2. Octave-XML, developed by Thomas Geiger, from Austria?.
Also people start talking about the GUI frontier for Octave, and
John Swensen explains how he runs Octave interpreter with the Octave GUI properly. I wonder if writing an evaluating parse tree is a solution.

Nice to see the ideas flowing, and what makes Octave the really the diverse community on a shared platform it is. Being, Indian, I can appreciate the micro-worlds in macro-world.

Cheers,
Muthu

Octave-2-Matlab & Back (2)

Doxygen is such a nice tool. I cannot over emphasize its utility when trying to read a large codebase like GNU Octave's. So in trying to solve our parse-tree walker to convert Octave(Matlab) code to the other Matlab(Octave), we use the inheritance class list for GNU Octave Interpreter.

The following, is the class inheritance for the nodes on the parse-tree, that are in the AST. Now the interpreter or a conversion tool like Octave-2-Matlab will walk this tree. The AST Walker implements the visitor pattern supported by the tree object.



We can now look into the organization of the AST Walker, which is neatly documented/illustrated in this working piece of code at /cvs/octave/src/pt-pr-code.{cc|h}. This AST walker inherits from the base class as follows,

with tree_walker, a pure virtual base class, that encourages polymorphic design. Some of the members of tree_walker include

Now a debugger, code conversion tool, lint-like checker or whatever has to conveniently invoke the parser for the AST, and then walk it for whatever operations it needs to perform on the tree.

virtual void visit_index_expression (tree_index_expression &)=0
virtual void visit_matrix (tree_matrix &)=0
virtual void visit_cell (tree_cell &)=0
virtual void visit_multi_assignment (tree_multi_assignment &)=0
virtual void visit_no_op_command (tree_no_op_command &)=0
virtual void visit_constant (tree_constant &)=0
virtual void visit_fcn_handle (tree_fcn_handle &)=0
virtual void visit_parameter_list (tree_parameter_list &)=0

and many more functions that correspond to each language-object that derives from tree, and relevant functions that 'execute' that command or expression.

So the organization of all this code would start at the front end like:

   1 main
2 begin
3 AST=call-parser()
4 cwalk=Code_Converting_Walker(OCT2MAT)
5 cwalk.set_tree(AST)
6 cwalk.visit()
7 cleanup abracadabra
8 end
So that settles it all, if you know what to write within the pure virtual functions, visit_index_expression (tree_index_expression &)=0 and such. This conversion code is put simply in the previous blogpost by walking this tree in a post-order fashion.

Next time we could see how to write the Oct2Mat converter based on the parse tree, or atleast drive the tree-printer class code.

Cheers,
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

Monday, June 25, 2007

Octave-2-Matlab and back

Octave-Lint (3)
A few more additional features that could be copied from existing tools in a similar areas. Immediately we can think of all the warnings that GCC would spit out, when you compile C/C++ code using 'gcc -Wall code.cc'. Specifically, flagging unused variables, data type conversions in known stdlibrary of Octave functions and such. Its time, I could show the code, instead of waxing eloquent about it. Just that I dont have anything to show.

Octave -> Matlab -> Octave
Next in a similar spirit of the AST walker static-checker we could also have a program transformation tool for converting Octave -> Matlab code. I could 'ask' for type inference whenever it is 'confused' by 'ambiguous' nature of the code. This is exciting, in terms of the possibilities it opens up for us. It is going to be a long project, but defnitely rewarding one.
The Octave parser is widely acknowledged on our mailing lists as a superset of Matlab
language. Using the Octave parser, you could in theory convert Matlab code the other-way
round too. Now this is double bonanza, a kind of buy-1 and you get-1 free. Am I excited?

The starting point of the program transformation tool would be from Paul Kienzle's work on
Oct2Mat conversion script. A few things are simple to understand, and require plain substituion. Paul's code does it using an AWK script, and clearly mentions not wanting to Octave's parse tree for doing the conversion. Paul attempts to do it using filters in AWK, and aims for an independent program for doing the conversion. Since our aims are to take advantage of whatever we have (GPL, Octave, Parser) then we can stand on the shoulders of giants instead of 'lets start at the very beginning, Do Re Mi'. Paul's code helps us easily identify the transformation rules to be applied.

The Octave->Matlab transformation rules, applied in Paul's code include
  1. gsub("#" , "%"); convert # to %
  2. gsub("[(][)]",""); convert () to ' ' as Matlab 4.0 and Octave-2.1.50 dont support empty arguments. This is not a problem anymore.
  3. gsub("gset[^;%]*;",""); a graphics command conversion.
  4. gsub("gset[^;]*%","%"); a graphics command conversion.
  5. gsub("gset[^;]*$",""); a graphics command conversion.
  6. gsub("endfunction",""); Octave has endif, endfor, endfunction etc, so those get converted to end, or in Matlab nothing at all. Infact Matlab complains if you have the end keyword at the end of a function file.
  7. gsub("endif","end"); Same as 6.
  8. gsub("endwhile","end"); Same as 6.
  9. gsub("endfor","end"); Same as 6.
  10. gsub("end_try_catch","end"); Same as 6.
  11. gsub(/&&/,"\\&"); Matlab didnot have the short circuit logical operators initially, so we had to use the logical &.
  12. gsub("SEEK_CUR",0); They also seemingly lacked options for STDIO processing.
  13. gsub("SEEK_END",1); Same as 12.
  14. gsub("SEEK_SET",-1); Same as 12.
  15. gsub("usage","error"); Our usage() and print_usage() are replaced by error() function.
  16. gsub("__error_text__","lasterr"); Some Octave specific code
  17. gsub("unwind_protect_cleanup",""); Again Octave was first here, w.r.t try/catch ...
  18. gsub("end_unwind_protect",""); ... doing it the LISP way.
  19. gsub("unwind_protect",""); Same as [17, 18].
  20. gsub(/\|\|/,"|"); Logical short-circuit OR operator was NOT there in Matlab first. We got here earlier.
  21. gsub("!","~"); They dont know whats Not !
  22. gsub("[\\]$","..."); Line continuations are OK for us. Not them!

These rules are easily loaded as actions for a AST walker, for operators in Octave. We just generate the do_convert_oct_to_mat() on each node. As an example, the node for the operator short-circuit && will have the conversion rule embedded in this routine as follows,

string short_ckt_and::do_convert_oct_to_mat() {
l_str=left_node.do_convert_oct_to_mat();
r_str=right_node.do_convert_oct_to_mat();
str=l_str + " & " + r_str;
return str;
}

So basically it would be a really syntactically correct way of doing things for bulilding transformation tools from the existing code base. As explained earlied since Octave being a superset of Matlab we could also turn the world the other-way round.

Ofcourse the validity of the converter can be verified by doing the rountripping, O-2-M-2-O and running the code to see if we have the regression tests validated. John Eaton (JWE), has a set of regression cases for parser, which is again something we can re-use.

It would be fun to see this happen, and I suspect this is much lower hanging fruit than writing an evaluating AST-walker for building the profiler or a debugger.

Cheers,
Muthu

Creative Commons License