Treelayout and CGS

classic Classic list List threaded Threaded
12 messages Options
Reply | Threaded
Open this post in threaded view
|

Treelayout and CGS

salac.r
This post was updated on .
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Treelayout and CGS

David Bateman-3
R! wrote:

> Hi,
>
> As a part of my studies at Brno University of Technology -Faculty of
> Information Technology I am supposed to
> write some piece of code. I chose to participate in Octave project under
> guidance of Ivana Varekova.
> I have completed procedure Treelayout and CGS for now and hope it will be
> fine for you.All suggestions are welcome.
>
> I am looking forward to hearing from you
>                                                     Radek Salac
>
> http://www.nabble.com/file/p20446040/cgs.m cgs.m
> http://www.nabble.com/file/p20446040/treelayout.m treelayout.m
>  
Hi Radek,

It seems you aren't subscribed to the maintainers mailing list, as your
e-mail dated the 11th Nov was just received on the mailing list. Perhaps
you might like to subscribe at

https://www-old.cae.wisc.edu/mailman/listinfo/octave-maintainers

As for the code, this is great and the only comments I have are more
style issues. I can update the style and commit these if you like, or
you might like to read the "Octave Sources" section of

http://hg.savannah.gnu.org/hgweb/octave/file/tip/doc/interpreter/contrib.txi

that discusses some of the stylistic constraints on the Octave sources..
Tell me what you want to do, and I'll go ahead and commit these additions.

One point however, the utility of the "cgs" function is limited by the
fact that Octave doesn't yet have good pre-conditioner functions. It
would be great if that point was also addressed.

Regards
David


--
David Bateman                                [hidden email]
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

Reply | Threaded
Open this post in threaded view
|

Re: Treelayout and CGS

David Bateman-3
David Bateman wrote:
>
> http://hg.savannah.gnu.org/hgweb/octave/file/tip/doc/interpreter/contrib.txi 
>

One point not address in that style guide, is that we prefer not to have
CamelCase variable names in the Octave sources. So variables like
ChildNumber in your code would preferably be written as child_number.

regards
David

--
David Bateman                                [hidden email]
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

Reply | Threaded
Open this post in threaded view
|

Re: Treelayout and CGS

salac.r
This post was updated on .
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Treelayout and CGS

Jaroslav Hajek-2
On Fri, Nov 14, 2008 at 12:59 PM, R! <[hidden email]> wrote:

>
> Hi
>
> Thanks for your suggestion, I will look at the style issue and change the
> source code as it is mentioned in the documentation (sorry that I didn't do
> it before). But as I think that there will be also another suggestion I will
> wait some time with a new version.
>
> Radek Salac
>

Hi Radek,

this looks very nice. I have two comments to cgs:

1. Since calling `error' enters an exceptional state, no more function
code is executed (unless it's inside unwind_protect or try/catch
block). Therefore, you can write the error checks as
if (check1)
  error ("error1")
elseif (check2)
  error ("error2")
endif

... normal code ...

this saves you one indentation of the normal code.

2. What is the reason for supplying M0 & M1, if they are multiplied
anyway? I don't remember much from Krylov methods, but I guess perhaps
the purpose of specifying 2 preconditioners is so that their
multiplication can be avoided?


--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
Reply | Threaded
Open this post in threaded view
|

Re: Treelayout and CGS

David Bateman-3
Jaroslav Hajek wrote:
> 2. What is the reason for supplying M0 & M1, if they are multiplied
> anyway? I don't remember much from Krylov methods, but I guess perhaps
> the purpose of specifying 2 preconditioners is so that their
> multiplication can be avoided?
>
>
>  
The book at

ftp://ftp.netlib.org/templates/templates.ps

might help for this

Regards
David

--
David Bateman                                [hidden email]
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

Reply | Threaded
Open this post in threaded view
|

Re: Treelayout and CGS

salac.r
This post was updated on .
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Treelayout and CGS

David Bateman-2
R! wrote:
> Hi,
> So there is a second version, i clean cgs procedure so now it fit more the
> codding standards. I also find one bug which is solved now.
> In treelayout I also made some cleanup but I don't wont change the variable
> name much.I know that it's bad to have CamelCase variable but these script
> is closely conected with treeplot where are also these variables.
> New version:
> http://www.nabble.com/file/p20591667/cgs.m cgs.m
> http://www.nabble.com/file/p20591667/treelayout.m treelayout.m


Why do you accept two preconditioners in cgs and then just multiply them
together? The reason to pass two preconditioners in this manner is that
if M1 and M2 are sparse then (M1 * M2) has a good chance of having a lot
of fill in and so M1 * (M2 * res) will therefore be much faster than (M1
*  M2) * res and use a lot less memory. If M1 and M2 are sparse then
please don't multiply them together in this fashion, but rather do
something like

   precon = [];
   if (nargin == 5)
     precon = M1;
   elseif (nargin > 5)
     if (isparse(M1) && issparse(M2)
       precon = @(x) M1 * (M2 * x);
     else
       precon = M1 * M2;
     endif
   endif

Apart from that I can't see any issues. If you are happy with the above
change I'll commit these when I have the chance.

Cheers
David


--
David Bateman                                [hidden email]
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)
Reply | Threaded
Open this post in threaded view
|

Re: Treelayout and CGS

salac.r
This post was updated on .
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Treelayout and CGS

David Bateman-3
R! wrote:

> Well thanks for your opinion, I really appreciate  it.
> If you think that you can improve my code in any way please do it,but if you
> don't have time I will try it myself. But it's your idea so I think that you
> have the right to implement it :-).
>
> Thanks for all.
>
> Radek Salac.
>
>  
>
>  
Ok, I committed the attached changeset credited to you.

D.


--
David Bateman                                [hidden email]
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


# HG changeset patch
# User Radek Salac <[hidden email]>
# Date 1227276183 -3600
# Node ID 5d8461010f118881bd040c02e8bbebd6ba81ee77
# Parent  67d5e0fd3c7dffbc4fba9ada0a0399c423a14a25
Add the cgs and treelayout functions

diff --git a/doc/interpreter/contributors.in b/doc/interpreter/contributors.in
--- a/doc/interpreter/contributors.in
+++ b/doc/interpreter/contributors.in
@@ -170,6 +170,7 @@
 Olli Saarela
 Toni Saarela
 Juhani Saastamoinen
+Radek Salac
 Ben Sapp
 Alois Schloegl
 Michel D. Schmid
diff --git a/scripts/ChangeLog b/scripts/ChangeLog
--- a/scripts/ChangeLog
+++ b/scripts/ChangeLog
@@ -1,3 +1,8 @@
+2008-11-21  Radek Salac  <[hidden email]>
+
+ * sparse/cgs.m, sparse/treelayout.m: New functions.
+ * sparse/Makefile.in (SOURCES): Add them here.
+
 2008-11-14  David Bateman  <[hidden email]>
 
  * plot/__go_draw_axes__.m (do_tics_1): Support the minorick properties
diff --git a/scripts/sparse/Makefile.in b/scripts/sparse/Makefile.in
--- a/scripts/sparse/Makefile.in
+++ b/scripts/sparse/Makefile.in
@@ -32,10 +32,10 @@
 INSTALL_PROGRAM = @INSTALL_PROGRAM@
 INSTALL_DATA = @INSTALL_DATA@
 
-SOURCES = colperm.m etreeplot.m gplot.m nonzeros.m normest.m \
+SOURCES = cgs.m colperm.m etreeplot.m gplot.m nonzeros.m normest.m \
   pcg.m pcr.m spalloc.m spaugment.m spconvert.m spdiags.m speye.m \
   spfun.m sphcat.m spones.m sprand.m sprandn.m sprandsym.m spstats.m \
-  spvcat.m spy.m treeplot.m
+  spvcat.m spy.m treelayout.m treeplot.m
 
 DISTFILES = $(addprefix $(srcdir)/, Makefile.in $(SOURCES))
 
diff --git a/scripts/sparse/cgs.m b/scripts/sparse/cgs.m
new file mode 100644
--- /dev/null
+++ b/scripts/sparse/cgs.m
@@ -0,0 +1,160 @@
+## Copyright (C) 2008 Radek Salac
+##
+## 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 3 of the License, 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, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} cgs (@var{A}, @var{b})
+## @deftypefnx {Function File} {} cgs (@var{A}, @var{b}, @var{tol}, @var{maxit}, @var{M1}, @var{M2}, @var{x0})
+## This procedure attempts to solve a system of linear equations A*x = b for x.
+## The @var{A} must be square, symmetric and positive definite real matrix N*N.
+## The @var{b} must be a one column vector with a length of N.
+## The @var{tol} specifies the tolerance of the method, default value is 1e-6.
+## The @var{maxit} specifies the maximum number of iteration, default value is MIN(20,N).
+## The @var{M1} specifies a preconditioner, can also be a function handler which returns M\X.
+## The @var{M2} combined with @var{M1} defines preconditioner as preconditioner=M1*M2.
+## The @var{x0} is initial guess, default value is zeros(N,1).
+##
+## @end deftypefn
+
+function [x, flag, relres, iter, resvec] = cgs (A, b, tol, maxit, M1, M2, x0)
+
+  if (nargin < 2 || nargin > 7 || nargout > 5)
+    print_usage ();
+  elseif (!isnumeric (A) || rows (A) != columns (A))
+    error ("cgs: first argument must be a n-by-n matrix");
+  elseif (!isvector (b))
+    error ("cgs: b must be a vector");
+  elseif (rows (A) != rows (b))
+    error ("cgs: first and second argument must have the same number of rows");
+  elseif (nargin > 2 && !isscalar (tol))
+    error ("cgs: tol must be a scalar");
+  elseif (nargin > 3 && !isscalar (maxit))
+    error ("cgs: maxit must be a scalar");
+  elseif (nargin > 4 && ismatrix (M1) && (rows (M1) != rows (A) || columns (M1) != columns (A)))
+    error ("cgs: M1 must have the same number of rows and columns as A");
+  elseif (nargin > 5 && (!ismatrix (M2) || rows (M2) != rows (A) || columns (M2) != columns (A)))
+    error ("cgs: M2 must have the same number of rows and columns as A");
+  elseif (nargin > 6 && !isvector (x0))
+    error ("cgs: x0 must be a vector");
+  elseif (nargin > 6 && rows (x0) != rows (b))
+    error ("cgs: xO must have the same number of rows as b");
+  endif
+
+  ## default toleration
+  if (nargin < 3)
+    tol = 1e-6;
+  endif
+
+  ## default maximum number of iteration
+  if (nargin < 4)
+    maxit = min (rows (b),20);
+  endif
+
+
+  ## left preconditioner
+  precon = [];
+  if (nargin == 5)
+    precon = M1;
+  elseif (nargin > 5)
+    if (isparse(M1) && issparse(M2))
+      precon = @(x) M1 * (M2 * x);
+    else
+      precon = M1 * M2;
+    endif
+  endif
+
+  ## precon can by also function
+  if (nargin > 4 && isnumeric (precon))
+
+    ## we can compute inverse preconditioner and use quicker algorithm
+    if (det (precon) != 0)
+     precon=inv (precon);
+    else
+     error ("cgs: preconditioner is ill conditioned");
+    endif
+
+    ## we must make test if preconditioner isn't ill conditioned
+    if (isinf (cond (precon)));
+      error ("cgs: preconditioner is ill conditioned");
+    endif
+  endif
+
+  ## specifies initial estimate x0
+  if (nargin < 7)
+    x = zeros (rows (b), 1);
+  else
+    x = x0;
+  endif
+
+  relres = b - A * x;
+  ## vector of the residual norms for each iteration
+  resvec = [norm(relres)];
+  ro = 0;
+  norm_b = norm (b);
+  ## default behaviour we don't reach tolerance tol within maxit iterations
+  flag = 1;
+  for iter = 1 : maxit
+
+    if (nargin > 4 && isnumeric (precon))
+      ## we have computed inverse matrix so we can use quick algorithm
+      z = precon * relres;
+    elseif (nargin > 4)
+      ## our preconditioner is a function
+      z = feval (precon, relres);
+    else
+      ## we don't use preconditioning
+      z = relres;
+    endif
+
+    ## cache
+    ro_old = ro;
+    ro = relres' * z;
+    if (iter == 1)
+      p = z;
+    else
+      beta = ro / ro_old;
+      p = z + beta * p;
+    endif
+    q = A * p; #cache
+    alpha = ro / (p' * q);
+    x = x + alpha * p;
+
+    relres = relres - alpha * q;
+    resvec = [resvec; norm(relres)];
+
+    relres_distance = resvec (end) / norm_b;
+    if (relres_distance <= tol)
+      ## we reach tolerance tol within maxit iterations
+      flag = 0;
+      break;
+    elseif (resvec (end) == resvec (end - 1) )
+      ## the method stagnates
+      flag = 3;
+      break;
+    endif
+  endfor;
+
+  relres = relres_distance;
+endfunction
+
+
+
+%!demo
+%! % Solve system of A*x=b
+%! A=[5 -1 3;-1 2 -2;3 -2 3]
+%! b=[7;-1;4]
+%! [a,b,c,d,e]=cgs(A,b)
diff --git a/scripts/sparse/treelayout.m b/scripts/sparse/treelayout.m
new file mode 100644
--- /dev/null
+++ b/scripts/sparse/treelayout.m
@@ -0,0 +1,208 @@
+## Copyright (C) 2008 Ivana Varekova & Radek Salac
+##
+## 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 3 of the License, 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, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} treelayout (@var{Tree})
+## @deftypefnx {Function File} {} treelayout (@var{Tree}, @var{Permutation})
+## treelayout lays out a tree or a forest. The first argument @var{Tree} is a vector of
+## predecessors, optional parameter @var{Permutation} is an optional postorder permutation.
+## The complexity of the algorithm is O(n) in
+## terms of time and memory requirements.
+## @seealso{etreeplot, gplot,treeplot}
+## @end deftypefn
+
+function [XCoordinate, YCoordinate, Height, s] = treelayout (Tree, Permutation)
+  if (nargin < 1 || nargin > 2 || nargout > 4)
+    print_usage ();
+  elseif (! isvector (Tree) || rows (Tree) != 1 || ! isnumeric (Tree)
+        ||  any (Tree > length (Tree)) || any (Tree < 0) )
+    error ("treelayout: the first input argument must be a vector of predecessors");
+  else
+    ## make it a row vector
+    Tree = Tree(:)';
+
+    ## the count of nodes of the graph
+    NodNumber = length (Tree);
+    ## the number of children
+    ChildNumber = zeros (1, NodNumber + 1);
+
+
+    ## checking vector of predecessors
+    for i = 1 : NodNumber
+      if (Tree (i) < i)
+ ## this part of graph was checked before
+        continue;
+      endif
+
+      ## Try to find cicle in this part of graph
+      ## we use modified Floyd's cycle-finding algorithm
+      tortoise = Tree (i);
+      hare = Tree (tortoise);
+
+      while (tortoise != hare)
+ ## we end after find a cicle or when we reach a checked part of graph
+
+        if (hare < i)
+          ## this part of graph was checked before
+          break
+        endif
+
+        tortoise = Tree (tortoise);
+ ## hare will move faster than tortoise so in cicle hare
+ ## must reach tortoise
+        hare = Tree (Tree (hare));
+
+      endwhile
+
+      if (tortoise == hare)
+ ## if hare reach tortoise we find cicle
+        error ("treelayout: vector of predecessors has bad format");
+      endif
+
+    endfor
+    ## vector of predecessors has right format
+
+    for i = 1:NodNumber
+      ## VecOfChild is helping vector which is used to speed up the
+      ## choose of descendant nodes
+
+      ChildNumber (Tree (i) + 1) = ChildNumber (Tree (i) + 1) + 1;
+    endfor
+
+    Pos = 1;
+    for i = 1 : NodNumber + 1
+      Start (i) = Pos;
+      Help (i) = Pos;
+      Pos += ChildNumber (i);
+      Stop (i) = Pos;
+    endfor
+
+    if (nargin == 1)
+      for i = 1 : NodNumber
+        VecOfChild (Help (Tree (i) + 1)) = i;  
+        Help (Tree (i) + 1) = Help (Tree (i) + 1) + 1;
+      endfor
+    else
+      VecOfChild = Permutation;
+    endif
+
+
+    ## the number of "parent" (actual) node (it's descendants will be
+    ## browse in the next iteration)
+    ParNumber = 0;
+
+    ## the x-coordinate of the left most descendant of "parent node"
+    ## this value is increased in each leaf
+    LeftMost = 0;
+
+    ## the level of "parent" node (root level is NodNumber)
+    Level = NodNumber;
+
+    ## NodNumber - Max is the height of this graph
+    Max = NodNumber;
+
+    ## main stack - each item consists of two numbers - the number of
+    ## node and the number it's of parent node on the top of stack
+    ## there is "parent node"
+    St = [-1, 0];
+
+    #number of vertices s in the top-level separator
+    s = 0;
+    # flag which says if we are in top level separator
+    topLevel = 1;
+    ## the top of the stack
+    while (ParNumber != -1)
+      if (Start(ParNumber + 1) < Stop(ParNumber + 1))
+        idx = VecOfChild (Start (ParNumber + 1) : Stop (ParNumber + 1) - 1);
+      else
+        idx = zeros (1, 0);
+      endif
+
+      ## add to idx the vector of parent descendants
+      St = [St ; [idx', ones(fliplr(size(idx))) * ParNumber]];
+
+      # we are in top level separator when we have one children
+      ## and the flag is 1
+      if (columns(idx) == 1 && topLevel ==1 )
+        s += 1;
+      else
+        # we arent in top level separator now
+        topLevel = 0;
+      endif
+      ## if there is not any descendant of "parent node":
+      if (St(end,2) != ParNumber)
+       LeftMost = LeftMost + 1;
+       XCoordinateR(ParNumber) = LeftMost;          
+       Max = min (Max, Level);
+       if ((length(St) > 1) && (find((shift(St,1)-St) == 0) >1)
+   && St(end,2) != St(end-1,2))
+  ## return to the nearest branching the position to return
+  ## position is the position on the stack, where should be
+          ## started further search (there are two nodes which has the
+          ## same parent node)
+
+          Position = (find ((shift (St(:, 2), 1) - St(:, 2)) == 0))(end)+1;
+          ParNumberVec = St(Position : end, 2);
+
+          ## the vector of removed nodes (the content of stack form
+          ## position to end)
+
+          Level = Level + length(ParNumberVec);
+
+  ## the level have to be decreased
+
+          XCoordinateR(ParNumberVec) = LeftMost;
+          St(Position:end, :) = [];
+        endif
+
+        ## remove the next node from "searched branch"
+
+        St(end, :) = [];
+ ## choose new "parent node"
+        ParNumber = St(end, 1);
+ ## if there is another branch start to search it
+ if (ParNumber != -1)
+          YCoordinate(ParNumber) = Level;
+          XCoordinateL(ParNumber) = LeftMost + 1;
+ endif
+      else
+
+        ## there were descendants of "parent nod" choose the last of
+        ## them and go on through it
+        Level--;
+        ParNumber = St(end, 1);
+        YCoordinate(ParNumber) = Level;    
+        XCoordinateL(ParNumber) = LeftMost+1;
+      endif
+    endwhile
+
+    ## calculate the x coordinates (the known values are the position
+    ## of most left and most right descendants)
+    XCoordinate = (XCoordinateL + XCoordinateR) / 2;
+
+    Height = NodNumber - Max - 1;
+  endif
+endfunction
+
+%!demo
+%! % Compute a simple tree layout
+%! [x,y,h,s]=treelayout([0 1 2 2])
+
+%!demo
+%! % Compute a simple tree layout with defined postorder permutation
+%! [x,y,h,s]=treelayout([0 1 2 2],[1 2 3 4])
Reply | Threaded
Open this post in threaded view
|

Re: Treelayout and CGS

salac.r
This post was updated on .
CONTENTS DELETED
The author has deleted this message.
Reply | Threaded
Open this post in threaded view
|

Re: Treelayout and CGS

David Bateman-2
R! wrote:
> Hi,
>
> Thanks all for your help and advice.Can I ask you what I shall do now for
> accepting these script to Octave please?
>
> Radek Salac

They already are.. That is what I meant by committed in the previous
e-mail. You are also credited as the author of the changeset in th
mercurial logs and added to the contributors file and so will be
credited in a 3.1.52 and after version of the manual

Regards
David


--
David Bateman                                [hidden email]
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)