Re: Treelayout and CGS

[David Bateman]

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

# HG changeset patch
# User Radek Salac <address@hidden>
# 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  <address@hidden>
+
+       * sparse/cgs.m, sparse/treelayout.m: New functions.
+       * sparse/Makefile.in (SOURCES): Add them here.
+
 2008-11-14  David Bateman  <address@hidden>
 
        * 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])