--- a/bindings/python/dionysus/__init__.py Tue Jun 05 14:50:58 2012 -0700
+++ b/bindings/python/dionysus/__init__.py Tue Jun 05 17:04:10 2012 -0700
@@ -3,6 +3,7 @@
from zigzag import *
from adaptor import *
import circular
+import viewer
def init_with_none(self, iter, data = None): # convenience: data defaults to None
self._cpp_init_(iter, data)
--- a/bindings/python/dionysus/circular.py Tue Jun 05 14:50:58 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,59 +0,0 @@
-def smooth(filtration, cocycle):
- from cvxopt import spmatrix, matrix
- from cvxopt.blas import copy
- from lsqr import lsqr
-
- coefficient = []
- coface_indices = []
- face_indices = []
- for i,s in enumerate(filtration):
- if s.dimension() > 2: continue
-
- c = 1
- for sb in s.boundary:
- j = filtration(sb)
- coefficient.append(c)
- coface_indices.append(i)
- face_indices.append(j)
- c *= -1
-
- # D is a coboundary matrix
- dimension = max(max(coface_indices), max(face_indices)) + 1
- D = spmatrix(coefficient, coface_indices, face_indices, (dimension, dimension))
-
- z = spmatrix([zz[0] for zz in cocycle],
- [zz[1] for zz in cocycle],
- [0 for zz in cocycle], (dimension, 1))
-
- v1 = D * z
- # print "D^2 is zero:", not bool(D*D)
- # print "D*z is zero:", not bool(v1)
- z = matrix(z)
-
- def Dfun(x,y,trans = 'N'):
- if trans == 'N':
- copy(D * x, y)
- elif trans == 'T':
- copy(D.T * x, y)
- else:
- assert False, "Unexpected trans parameter"
-
- tol = 1e-10
- show = False
- maxit = None
- solution = lsqr(Dfun, matrix(z), show = show, atol = tol, btol = tol, itnlim = maxit)
-
- z_smooth = z - D*solution[0]
-
- # print sum(z_smooth**2)
- # assert sum((D*z_smooth)**2) < tol and sum((D.T*z_smooth)**2) < tol, "Expected a harmonic cocycle"
- if not (sum((D*z_smooth)**2) < tol and sum((D.T*z_smooth)**2) < tol):
- print "Expected a harmonic cocycle:", sum((D*z_smooth)**2), sum((D.T*z_smooth)**2)
-
- values = {}
- vertices = ((i,s) for (i,s) in enumerate(filtration) if s.dimension() == 0)
- for i,s in vertices:
- v = [v for v in s.vertices][0]
- values[v] = solution[0][i]
-
- return values
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/dionysus/circular/__init__.py Tue Jun 05 17:04:10 2012 -0700
@@ -0,0 +1,59 @@
+def smooth(filtration, cocycle):
+ from cvxopt import spmatrix, matrix
+ from cvxopt.blas import copy
+ from lsqr import lsqr
+
+ coefficient = []
+ coface_indices = []
+ face_indices = []
+ for i,s in enumerate(filtration):
+ if s.dimension() > 2: continue
+
+ c = 1
+ for sb in s.boundary:
+ j = filtration(sb)
+ coefficient.append(c)
+ coface_indices.append(i)
+ face_indices.append(j)
+ c *= -1
+
+ # D is a coboundary matrix
+ dimension = max(max(coface_indices), max(face_indices)) + 1
+ D = spmatrix(coefficient, coface_indices, face_indices, (dimension, dimension))
+
+ z = spmatrix([zz[0] for zz in cocycle],
+ [zz[1] for zz in cocycle],
+ [0 for zz in cocycle], (dimension, 1))
+
+ v1 = D * z
+ # print "D^2 is zero:", not bool(D*D)
+ # print "D*z is zero:", not bool(v1)
+ z = matrix(z)
+
+ def Dfun(x,y,trans = 'N'):
+ if trans == 'N':
+ copy(D * x, y)
+ elif trans == 'T':
+ copy(D.T * x, y)
+ else:
+ assert False, "Unexpected trans parameter"
+
+ tol = 1e-10
+ show = False
+ maxit = None
+ solution = lsqr(Dfun, matrix(z), show = show, atol = tol, btol = tol, itnlim = maxit)
+
+ z_smooth = z - D*solution[0]
+
+ # print sum(z_smooth**2)
+ # assert sum((D*z_smooth)**2) < tol and sum((D.T*z_smooth)**2) < tol, "Expected a harmonic cocycle"
+ if not (sum((D*z_smooth)**2) < tol and sum((D.T*z_smooth)**2) < tol):
+ print "Expected a harmonic cocycle:", sum((D*z_smooth)**2), sum((D.T*z_smooth)**2)
+
+ values = {}
+ vertices = ((i,s) for (i,s) in enumerate(filtration) if s.dimension() == 0)
+ for i,s in vertices:
+ v = [v for v in s.vertices][0]
+ values[v] = solution[0][i]
+
+ return values
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/dionysus/circular/lsqr.py Tue Jun 05 17:04:10 2012 -0700
@@ -0,0 +1,408 @@
+# LSQR solver from http://pages.cs.wisc.edu/~kline/cvxopt/
+
+from cvxopt import matrix
+from cvxopt.lapack import *
+from cvxopt.blas import *
+from math import sqrt
+
+"""
+a,b are scalars
+
+On exit, returns scalars c,s,r
+"""
+def SymOrtho(a,b):
+ aa=abs(a)
+ ab=abs(b)
+ if b==0.:
+ s=0.
+ r=aa
+ if aa==0.:
+ c=1.
+ else:
+ c=a/aa
+ elif a==0.:
+ c=0.
+ s=b/ab
+ r=ab
+ elif ab>=aa:
+ sb=1
+ if b<0: sb=-1
+ tau=a/b
+ s=sb*(1+tau**2)**-0.5
+ c=s*tau
+ r=b/s
+ elif aa>ab:
+ sa=1
+ if a<0: sa=-1
+ tau=b/a
+ c=sa*(1+tau**2)**-0.5
+ s=c*tau
+ r=a/c
+
+ return c,s,r
+
+"""
+
+It is usually recommended to use SYMMLQ for symmetric matrices
+
+Requires the syntax
+ A(x,y) == y:=[A]*x
+and
+ A(x,y,trans='T') == y:=[A.T]*x
+
+comments with '###' are followed by the intent of the original matlab
+code. This may be useful for debugging.
+
+"""
+
+def lsqr( A, b, damp=0.0, atol=1e-8, btol=1e-8, conlim=1e8, itnlim=None, show=False, wantvar=False):
+ """
+
+ [ x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var ]...
+ = lsqr( m, n, 'aprod', iw, rw, b, damp, atol, btol, conlim, itnlim, show );
+
+ LSQR solves Ax = b or min ||b - Ax||_2 if damp = 0,
+ or min || (b) - ( A )x || otherwise.
+ || (0) (damp I) ||2
+ A is an m by n matrix defined by y = aprod( mode,m,n,x,iw,rw ),
+ where the parameter 'aprodname' refers to a function 'aprod' that
+ performs the matrix-vector operations.
+ If mode = 1, aprod must return y = Ax without altering x.
+ If mode = 2, aprod must return y = A'x without altering x.
+ WARNING: The file containing the function 'aprod'
+ must not be called aprodname.m !!!!
+
+ -----------------------------------------------------------------------
+ LSQR uses an iterative (conjugate-gradient-like) method.
+ For further information, see
+ 1. C. C. Paige and M. A. Saunders (1982a).
+ LSQR: An algorithm for sparse linear equations and sparse least squares,
+ ACM TOMS 8(1), 43-71.
+ 2. C. C. Paige and M. A. Saunders (1982b).
+ Algorithm 583. LSQR: Sparse linear equations and least squares problems,
+ ACM TOMS 8(2), 195-209.
+ 3. M. A. Saunders (1995). Solution of sparse rectangular systems using
+ LSQR and CRAIG, BIT 35, 588-604.
+
+ Input parameters:
+ iw, rw are not used by lsqr, but are passed to aprod.
+ atol, btol are stopping tolerances. If both are 1.0e-9 (say),
+ the final residual norm should be accurate to about 9 digits.
+ (The final x will usually have fewer correct digits,
+ depending on cond(A) and the size of damp.)
+ conlim is also a stopping tolerance. lsqr terminates if an estimate
+ of cond(A) exceeds conlim. For compatible systems Ax = b,
+ conlim could be as large as 1.0e+12 (say). For least-squares
+ problems, conlim should be less than 1.0e+8.
+ Maximum precision can be obtained by setting
+ atol = btol = conlim = zero, but the number of iterations
+ may then be excessive.
+ itnlim is an explicit limit on iterations (for safety).
+ show = 1 gives an iteration log,
+ show = 0 suppresses output.
+
+ Output parameters:
+ x is the final solution.
+ istop gives the reason for termination.
+ istop = 1 means x is an approximate solution to Ax = b.
+ = 2 means x approximately solves the least-squares problem.
+ r1norm = norm(r), where r = b - Ax.
+ r2norm = sqrt( norm(r)^2 + damp^2 * norm(x)^2 )
+ = r1norm if damp = 0.
+ anorm = estimate of Frobenius norm of Abar = [ A ].
+ [damp*I]
+ acond = estimate of cond(Abar).
+ arnorm = estimate of norm(A'*r - damp^2*x).
+ xnorm = norm(x).
+ var (if present) estimates all diagonals of (A'A)^{-1} (if damp=0)
+ or more generally (A'A + damp^2*I)^{-1}.
+ This is well defined if A has full column rank or damp > 0.
+ (Not sure what var means if rank(A) < n and damp = 0.)
+
+
+ 1990: Derived from Fortran 77 version of LSQR.
+ 22 May 1992: bbnorm was used incorrectly. Replaced by anorm.
+ 26 Oct 1992: More input and output parameters added.
+ 01 Sep 1994: Matrix-vector routine is now a parameter 'aprodname'.
+ Print log reformatted.
+ 14 Jun 1997: show added to allow printing or not.
+ 30 Jun 1997: var added as an optional output parameter.
+ 07 Aug 2002: Output parameter rnorm replaced by r1norm and r2norm.
+ Michael Saunders, Systems Optimization Laboratory,
+ Dept of MS&E, Stanford University.
+ -----------------------------------------------------------------------
+ """
+ """
+ Initialize.
+ """
+ n=len(b)
+ m=n
+ if itnlim is None: itnlim=2*n
+
+ msg=('The exact solution is x = 0 ',
+ 'Ax - b is small enough, given atol, btol ',
+ 'The least-squares solution is good enough, given atol ',
+ 'The estimate of cond(Abar) has exceeded conlim ',
+ 'Ax - b is small enough for this machine ',
+ 'The least-squares solution is good enough for this machine',
+ 'Cond(Abar) seems to be too large for this machine ',
+ 'The iteration limit has been reached ');
+
+ var = matrix(0.,(n,1));
+
+ if show:
+ print ' '
+ print 'LSQR Least-squares solution of Ax = b'
+ str1 = 'The matrix A has %8g rows and %8g cols' % (m, n)
+ str2 = 'damp = %20.14e wantvar = %8g' %( damp,wantvar)
+ str3 = 'atol = %8.2e conlim = %8.2e'%( atol, conlim)
+ str4 = 'btol = %8.2e itnlim = %8g' %( btol, itnlim)
+ print str1
+ print str2
+ print str3
+ print str4
+
+ itn = 0; istop = 0; nstop = 0;
+ ctol = 0;
+ if conlim > 0: ctol = 1/conlim
+ anorm = 0; acond = 0;
+ dampsq = damp**2; ddnorm = 0; res2 = 0;
+ xnorm = 0; xxnorm = 0; z = 0;
+ cs2 = -1; sn2 = 0;
+
+ """
+ Set up the first vectors u and v for the bidiagonalization.
+ These satisfy beta*u = b, alfa*v = A'u.
+ """
+ __x = matrix(0., (n,1)) # a matrix for temporary holding
+ v = matrix(0., (n,1))
+ u = +b;
+ x = matrix(0., (n,1))
+ alfa = 0;
+ beta = nrm2( u );
+ w = matrix(0., (n,1))
+
+ if beta > 0:
+ ### u = (1/beta) * u;
+ ### v = feval( aprodname, 2, m, n, u, iw, rw );
+ scal(1/beta,u)
+ A(u,v,trans='T'); #v = feval( aprodname, 2, m, n, u, iw, rw );
+ alfa = nrm2( v );
+
+ if alfa > 0:
+ ### v = (1/alfa) * v;
+ scal(1/alfa,v)
+ copy(v,w)
+
+
+ rhobar = alfa; phibar = beta; bnorm = beta;
+ rnorm = beta;
+ r1norm = rnorm;
+ r2norm = rnorm;
+
+ # reverse the order here from the original matlab code because
+ # there was an error on return when arnorm==0
+ arnorm = alfa * beta;
+ if arnorm == 0:
+ print msg[0];
+ return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var
+
+ head1 = ' Itn x[0] r1norm r2norm ';
+ head2 = ' Compatible LS Norm A Cond A';
+
+ if show:
+ print ' '
+ print head1, head2
+ test1 = 1; test2 = alfa / beta;
+ str1 = '%6g %12.5e' %( itn, x[0] );
+ str2 = ' %10.3e %10.3e'%( r1norm, r2norm );
+ str3 = ' %8.1e %8.1e' %( test1, test2 );
+ print str1, str2, str3
+
+ """
+ %------------------------------------------------------------------
+ % Main iteration loop.
+ %------------------------------------------------------------------
+ """
+ while itn < itnlim:
+ itn = itn + 1;
+ """
+ % Perform the next step of the bidiagonalization to obtain the
+ % next beta, u, alfa, v. These satisfy the relations
+ % beta*u = a*v - alfa*u,
+ % alfa*v = A'*u - beta*v.
+ """
+ ### u = feval( aprodname, 1, m, n, v, iw, rw ) - alfa*u;
+ copy(u, __x)
+ A(v,u)
+ axpy(__x,u,-alfa)
+
+ beta = nrm2( u );
+ if beta > 0:
+ ### u = (1/beta) * u;
+ scal(1/beta,u)
+ anorm = sqrt(anorm**2 + alfa**2 + beta**2 + damp**2);
+ ### v = feval( aprodname, 2, m, n, u, iw, rw ) - beta*v;
+ copy(v,__x)
+ A(u,v,trans='T')
+ axpy(__x,v,-beta)
+
+ alfa = nrm2( v );
+ if alfa > 0:
+ ### v = (1/alfa) * v;
+ scal(1/alfa, v)
+
+ """
+ % Use a plane rotation to eliminate the damping parameter.
+ % This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
+ """
+
+ rhobar1 = sqrt(rhobar**2 + damp**2);
+ cs1 = rhobar / rhobar1;
+ sn1 = damp / rhobar1;
+ psi = sn1 * phibar;
+ phibar = cs1 * phibar;
+ """
+ % Use a plane rotation to eliminate the subdiagonal element (beta)
+ % of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
+ """
+
+
+ ###cs = rhobar1/ rho;
+ ###sn = beta / rho;
+ cs,sn,rho = SymOrtho(rhobar1,beta)
+
+ theta = sn * alfa;
+ rhobar = - cs * alfa;
+ phi = cs * phibar;
+ phibar = sn * phibar;
+ tau = sn * phi;
+ """
+ % Update x and w.
+ """
+ t1 = phi /rho;
+ t2 = - theta/rho;
+ dk = (1/rho)*w;
+
+ ### x = x + t1*w;
+ axpy(w,x,t1)
+ ### w = v + t2*w;
+ scal(t2,w)
+ axpy(v,w)
+ ddnorm = ddnorm + nrm2(dk)**2;
+ if wantvar:
+ ### var = var + dk.*dk;
+ axpy(dk**2, var)
+ """
+ % Use a plane rotation on the right to eliminate the
+ % super-diagonal element (theta) of the upper-bidiagonal matrix.
+ % Then use the result to estimate norm(x).
+ """
+
+ delta = sn2 * rho;
+ gambar = - cs2 * rho;
+ rhs = phi - delta * z;
+ zbar = rhs / gambar;
+ xnorm = sqrt(xxnorm + zbar**2);
+ gamma = sqrt(gambar**2 +theta**2);
+ cs2 = gambar / gamma;
+ sn2 = theta / gamma;
+ z = rhs / gamma;
+ xxnorm = xxnorm + z**2;
+ """
+ % Test for convergence.
+ % First, estimate the condition of the matrix Abar,
+ % and the norms of rbar and Abar'rbar.
+ """
+ acond = anorm * sqrt(ddnorm);
+ res1 = phibar**2;
+ res2 = res2 + psi**2;
+ rnorm = sqrt( res1 + res2 );
+ arnorm = alfa * abs( tau );
+ """
+ % 07 Aug 2002:
+ % Distinguish between
+ % r1norm = ||b - Ax|| and
+ % r2norm = rnorm in current code
+ % = sqrt(r1norm^2 + damp^2*||x||^2).
+ % Estimate r1norm from
+ % r1norm = sqrt(r2norm^2 - damp^2*||x||^2).
+ % Although there is cancellation, it might be accurate enough.
+ """
+ r1sq = rnorm**2 - dampsq * xxnorm;
+ r1norm = sqrt( abs(r1sq) );
+ if r1sq < 0: r1norm = - r1norm;
+ r2norm = rnorm;
+ """
+ % Now use these norms to estimate certain other quantities,
+ % some of which will be small near a solution.
+ """
+ test1 = rnorm / bnorm;
+ test2 = arnorm/( anorm * rnorm );
+ test3 = 1 / acond;
+ t1 = test1 / (1 + anorm * xnorm / bnorm);
+ rtol = btol + atol * anorm * xnorm / bnorm;
+ """
+ % The following tests guard against extremely small values of
+ % atol, btol or ctol. (The user may have set any or all of
+ % the parameters atol, btol, conlim to 0.)
+ % The effect is equivalent to the normal tests using
+ % atol = eps, btol = eps, conlim = 1/eps.
+ """
+ if itn >= itnlim : istop = 7;
+ if 1 + test3 <= 1: istop = 6;
+ if 1 + test2 <= 1: istop = 5;
+ if 1 + t1 <= 1: istop = 4;
+ """
+ % Allow for tolerances set by the user.
+ """
+ if test3 <= ctol: istop = 3;
+ if test2 <= atol: istop = 2;
+ if test1 <= rtol: istop = 1;
+ """
+ % See if it is time to print something.
+ """
+ prnt = False;
+ if n <= 40 : prnt = True;
+ if itn <= 10 : prnt = True;
+ if itn >= itnlim-10: prnt = True;
+ # if itn%10 == 0 : prnt = True;
+ if test3 <= 2*ctol : prnt = True;
+ if test2 <= 10*atol : prnt = True;
+ if test1 <= 10*rtol : prnt = True;
+ if istop != 0 : prnt = True;
+
+ if prnt:
+ if show:
+ str1 = '%6g %12.5e'% (itn, x[0] );
+ str2 = ' %10.3e %10.3e'% (r1norm, r2norm );
+ str3 = ' %8.1e %8.1e'% ( test1, test2 );
+ str4 = ' %8.1e %8.1e'% ( anorm, acond );
+ print str1, str2, str3, str4
+
+ if istop != 0: break
+
+ """
+ % End of iteration loop.
+ % Print the stopping condition.
+ """
+ if show:
+ print ' '
+ print 'LSQR finished'
+ print msg[istop]
+ print ' '
+ str1 = 'istop =%8g r1norm =%8.1e'% ( istop, r1norm );
+ str2 = 'anorm =%8.1e arnorm =%8.1e'%( anorm, arnorm );
+ str3 = 'itn =%8g r2norm =%8.1e'% ( itn, r2norm );
+ str4 = 'acond =%8.1e xnorm =%8.1e'%( acond, xnorm );
+ print str1+ ' ' +str2
+ print str3+ ' ' +str4
+ print ' '
+
+ return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var
+
+ """
+ %-----------------------------------------------------------------------
+ % End of lsqr.m
+ %-----------------------------------------------------------------------
+ """
--- a/bindings/python/dionysus/lsqr.py Tue Jun 05 14:50:58 2012 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,408 +0,0 @@
-# LSQR solver from http://pages.cs.wisc.edu/~kline/cvxopt/
-
-from cvxopt import matrix
-from cvxopt.lapack import *
-from cvxopt.blas import *
-from math import sqrt
-
-"""
-a,b are scalars
-
-On exit, returns scalars c,s,r
-"""
-def SymOrtho(a,b):
- aa=abs(a)
- ab=abs(b)
- if b==0.:
- s=0.
- r=aa
- if aa==0.:
- c=1.
- else:
- c=a/aa
- elif a==0.:
- c=0.
- s=b/ab
- r=ab
- elif ab>=aa:
- sb=1
- if b<0: sb=-1
- tau=a/b
- s=sb*(1+tau**2)**-0.5
- c=s*tau
- r=b/s
- elif aa>ab:
- sa=1
- if a<0: sa=-1
- tau=b/a
- c=sa*(1+tau**2)**-0.5
- s=c*tau
- r=a/c
-
- return c,s,r
-
-"""
-
-It is usually recommended to use SYMMLQ for symmetric matrices
-
-Requires the syntax
- A(x,y) == y:=[A]*x
-and
- A(x,y,trans='T') == y:=[A.T]*x
-
-comments with '###' are followed by the intent of the original matlab
-code. This may be useful for debugging.
-
-"""
-
-def lsqr( A, b, damp=0.0, atol=1e-8, btol=1e-8, conlim=1e8, itnlim=None, show=False, wantvar=False):
- """
-
- [ x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var ]...
- = lsqr( m, n, 'aprod', iw, rw, b, damp, atol, btol, conlim, itnlim, show );
-
- LSQR solves Ax = b or min ||b - Ax||_2 if damp = 0,
- or min || (b) - ( A )x || otherwise.
- || (0) (damp I) ||2
- A is an m by n matrix defined by y = aprod( mode,m,n,x,iw,rw ),
- where the parameter 'aprodname' refers to a function 'aprod' that
- performs the matrix-vector operations.
- If mode = 1, aprod must return y = Ax without altering x.
- If mode = 2, aprod must return y = A'x without altering x.
- WARNING: The file containing the function 'aprod'
- must not be called aprodname.m !!!!
-
- -----------------------------------------------------------------------
- LSQR uses an iterative (conjugate-gradient-like) method.
- For further information, see
- 1. C. C. Paige and M. A. Saunders (1982a).
- LSQR: An algorithm for sparse linear equations and sparse least squares,
- ACM TOMS 8(1), 43-71.
- 2. C. C. Paige and M. A. Saunders (1982b).
- Algorithm 583. LSQR: Sparse linear equations and least squares problems,
- ACM TOMS 8(2), 195-209.
- 3. M. A. Saunders (1995). Solution of sparse rectangular systems using
- LSQR and CRAIG, BIT 35, 588-604.
-
- Input parameters:
- iw, rw are not used by lsqr, but are passed to aprod.
- atol, btol are stopping tolerances. If both are 1.0e-9 (say),
- the final residual norm should be accurate to about 9 digits.
- (The final x will usually have fewer correct digits,
- depending on cond(A) and the size of damp.)
- conlim is also a stopping tolerance. lsqr terminates if an estimate
- of cond(A) exceeds conlim. For compatible systems Ax = b,
- conlim could be as large as 1.0e+12 (say). For least-squares
- problems, conlim should be less than 1.0e+8.
- Maximum precision can be obtained by setting
- atol = btol = conlim = zero, but the number of iterations
- may then be excessive.
- itnlim is an explicit limit on iterations (for safety).
- show = 1 gives an iteration log,
- show = 0 suppresses output.
-
- Output parameters:
- x is the final solution.
- istop gives the reason for termination.
- istop = 1 means x is an approximate solution to Ax = b.
- = 2 means x approximately solves the least-squares problem.
- r1norm = norm(r), where r = b - Ax.
- r2norm = sqrt( norm(r)^2 + damp^2 * norm(x)^2 )
- = r1norm if damp = 0.
- anorm = estimate of Frobenius norm of Abar = [ A ].
- [damp*I]
- acond = estimate of cond(Abar).
- arnorm = estimate of norm(A'*r - damp^2*x).
- xnorm = norm(x).
- var (if present) estimates all diagonals of (A'A)^{-1} (if damp=0)
- or more generally (A'A + damp^2*I)^{-1}.
- This is well defined if A has full column rank or damp > 0.
- (Not sure what var means if rank(A) < n and damp = 0.)
-
-
- 1990: Derived from Fortran 77 version of LSQR.
- 22 May 1992: bbnorm was used incorrectly. Replaced by anorm.
- 26 Oct 1992: More input and output parameters added.
- 01 Sep 1994: Matrix-vector routine is now a parameter 'aprodname'.
- Print log reformatted.
- 14 Jun 1997: show added to allow printing or not.
- 30 Jun 1997: var added as an optional output parameter.
- 07 Aug 2002: Output parameter rnorm replaced by r1norm and r2norm.
- Michael Saunders, Systems Optimization Laboratory,
- Dept of MS&E, Stanford University.
- -----------------------------------------------------------------------
- """
- """
- Initialize.
- """
- n=len(b)
- m=n
- if itnlim is None: itnlim=2*n
-
- msg=('The exact solution is x = 0 ',
- 'Ax - b is small enough, given atol, btol ',
- 'The least-squares solution is good enough, given atol ',
- 'The estimate of cond(Abar) has exceeded conlim ',
- 'Ax - b is small enough for this machine ',
- 'The least-squares solution is good enough for this machine',
- 'Cond(Abar) seems to be too large for this machine ',
- 'The iteration limit has been reached ');
-
- var = matrix(0.,(n,1));
-
- if show:
- print ' '
- print 'LSQR Least-squares solution of Ax = b'
- str1 = 'The matrix A has %8g rows and %8g cols' % (m, n)
- str2 = 'damp = %20.14e wantvar = %8g' %( damp,wantvar)
- str3 = 'atol = %8.2e conlim = %8.2e'%( atol, conlim)
- str4 = 'btol = %8.2e itnlim = %8g' %( btol, itnlim)
- print str1
- print str2
- print str3
- print str4
-
- itn = 0; istop = 0; nstop = 0;
- ctol = 0;
- if conlim > 0: ctol = 1/conlim
- anorm = 0; acond = 0;
- dampsq = damp**2; ddnorm = 0; res2 = 0;
- xnorm = 0; xxnorm = 0; z = 0;
- cs2 = -1; sn2 = 0;
-
- """
- Set up the first vectors u and v for the bidiagonalization.
- These satisfy beta*u = b, alfa*v = A'u.
- """
- __x = matrix(0., (n,1)) # a matrix for temporary holding
- v = matrix(0., (n,1))
- u = +b;
- x = matrix(0., (n,1))
- alfa = 0;
- beta = nrm2( u );
- w = matrix(0., (n,1))
-
- if beta > 0:
- ### u = (1/beta) * u;
- ### v = feval( aprodname, 2, m, n, u, iw, rw );
- scal(1/beta,u)
- A(u,v,trans='T'); #v = feval( aprodname, 2, m, n, u, iw, rw );
- alfa = nrm2( v );
-
- if alfa > 0:
- ### v = (1/alfa) * v;
- scal(1/alfa,v)
- copy(v,w)
-
-
- rhobar = alfa; phibar = beta; bnorm = beta;
- rnorm = beta;
- r1norm = rnorm;
- r2norm = rnorm;
-
- # reverse the order here from the original matlab code because
- # there was an error on return when arnorm==0
- arnorm = alfa * beta;
- if arnorm == 0:
- print msg[0];
- return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var
-
- head1 = ' Itn x[0] r1norm r2norm ';
- head2 = ' Compatible LS Norm A Cond A';
-
- if show:
- print ' '
- print head1, head2
- test1 = 1; test2 = alfa / beta;
- str1 = '%6g %12.5e' %( itn, x[0] );
- str2 = ' %10.3e %10.3e'%( r1norm, r2norm );
- str3 = ' %8.1e %8.1e' %( test1, test2 );
- print str1, str2, str3
-
- """
- %------------------------------------------------------------------
- % Main iteration loop.
- %------------------------------------------------------------------
- """
- while itn < itnlim:
- itn = itn + 1;
- """
- % Perform the next step of the bidiagonalization to obtain the
- % next beta, u, alfa, v. These satisfy the relations
- % beta*u = a*v - alfa*u,
- % alfa*v = A'*u - beta*v.
- """
- ### u = feval( aprodname, 1, m, n, v, iw, rw ) - alfa*u;
- copy(u, __x)
- A(v,u)
- axpy(__x,u,-alfa)
-
- beta = nrm2( u );
- if beta > 0:
- ### u = (1/beta) * u;
- scal(1/beta,u)
- anorm = sqrt(anorm**2 + alfa**2 + beta**2 + damp**2);
- ### v = feval( aprodname, 2, m, n, u, iw, rw ) - beta*v;
- copy(v,__x)
- A(u,v,trans='T')
- axpy(__x,v,-beta)
-
- alfa = nrm2( v );
- if alfa > 0:
- ### v = (1/alfa) * v;
- scal(1/alfa, v)
-
- """
- % Use a plane rotation to eliminate the damping parameter.
- % This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
- """
-
- rhobar1 = sqrt(rhobar**2 + damp**2);
- cs1 = rhobar / rhobar1;
- sn1 = damp / rhobar1;
- psi = sn1 * phibar;
- phibar = cs1 * phibar;
- """
- % Use a plane rotation to eliminate the subdiagonal element (beta)
- % of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
- """
-
-
- ###cs = rhobar1/ rho;
- ###sn = beta / rho;
- cs,sn,rho = SymOrtho(rhobar1,beta)
-
- theta = sn * alfa;
- rhobar = - cs * alfa;
- phi = cs * phibar;
- phibar = sn * phibar;
- tau = sn * phi;
- """
- % Update x and w.
- """
- t1 = phi /rho;
- t2 = - theta/rho;
- dk = (1/rho)*w;
-
- ### x = x + t1*w;
- axpy(w,x,t1)
- ### w = v + t2*w;
- scal(t2,w)
- axpy(v,w)
- ddnorm = ddnorm + nrm2(dk)**2;
- if wantvar:
- ### var = var + dk.*dk;
- axpy(dk**2, var)
- """
- % Use a plane rotation on the right to eliminate the
- % super-diagonal element (theta) of the upper-bidiagonal matrix.
- % Then use the result to estimate norm(x).
- """
-
- delta = sn2 * rho;
- gambar = - cs2 * rho;
- rhs = phi - delta * z;
- zbar = rhs / gambar;
- xnorm = sqrt(xxnorm + zbar**2);
- gamma = sqrt(gambar**2 +theta**2);
- cs2 = gambar / gamma;
- sn2 = theta / gamma;
- z = rhs / gamma;
- xxnorm = xxnorm + z**2;
- """
- % Test for convergence.
- % First, estimate the condition of the matrix Abar,
- % and the norms of rbar and Abar'rbar.
- """
- acond = anorm * sqrt(ddnorm);
- res1 = phibar**2;
- res2 = res2 + psi**2;
- rnorm = sqrt( res1 + res2 );
- arnorm = alfa * abs( tau );
- """
- % 07 Aug 2002:
- % Distinguish between
- % r1norm = ||b - Ax|| and
- % r2norm = rnorm in current code
- % = sqrt(r1norm^2 + damp^2*||x||^2).
- % Estimate r1norm from
- % r1norm = sqrt(r2norm^2 - damp^2*||x||^2).
- % Although there is cancellation, it might be accurate enough.
- """
- r1sq = rnorm**2 - dampsq * xxnorm;
- r1norm = sqrt( abs(r1sq) );
- if r1sq < 0: r1norm = - r1norm;
- r2norm = rnorm;
- """
- % Now use these norms to estimate certain other quantities,
- % some of which will be small near a solution.
- """
- test1 = rnorm / bnorm;
- test2 = arnorm/( anorm * rnorm );
- test3 = 1 / acond;
- t1 = test1 / (1 + anorm * xnorm / bnorm);
- rtol = btol + atol * anorm * xnorm / bnorm;
- """
- % The following tests guard against extremely small values of
- % atol, btol or ctol. (The user may have set any or all of
- % the parameters atol, btol, conlim to 0.)
- % The effect is equivalent to the normal tests using
- % atol = eps, btol = eps, conlim = 1/eps.
- """
- if itn >= itnlim : istop = 7;
- if 1 + test3 <= 1: istop = 6;
- if 1 + test2 <= 1: istop = 5;
- if 1 + t1 <= 1: istop = 4;
- """
- % Allow for tolerances set by the user.
- """
- if test3 <= ctol: istop = 3;
- if test2 <= atol: istop = 2;
- if test1 <= rtol: istop = 1;
- """
- % See if it is time to print something.
- """
- prnt = False;
- if n <= 40 : prnt = True;
- if itn <= 10 : prnt = True;
- if itn >= itnlim-10: prnt = True;
- # if itn%10 == 0 : prnt = True;
- if test3 <= 2*ctol : prnt = True;
- if test2 <= 10*atol : prnt = True;
- if test1 <= 10*rtol : prnt = True;
- if istop != 0 : prnt = True;
-
- if prnt:
- if show:
- str1 = '%6g %12.5e'% (itn, x[0] );
- str2 = ' %10.3e %10.3e'% (r1norm, r2norm );
- str3 = ' %8.1e %8.1e'% ( test1, test2 );
- str4 = ' %8.1e %8.1e'% ( anorm, acond );
- print str1, str2, str3, str4
-
- if istop != 0: break
-
- """
- % End of iteration loop.
- % Print the stopping condition.
- """
- if show:
- print ' '
- print 'LSQR finished'
- print msg[istop]
- print ' '
- str1 = 'istop =%8g r1norm =%8.1e'% ( istop, r1norm );
- str2 = 'anorm =%8.1e arnorm =%8.1e'%( anorm, arnorm );
- str3 = 'itn =%8g r2norm =%8.1e'% ( itn, r2norm );
- str4 = 'acond =%8.1e xnorm =%8.1e'%( acond, xnorm );
- print str1+ ' ' +str2
- print str3+ ' ' +str4
- print ' '
-
- return x, istop, itn, r1norm, r2norm, anorm, acond, arnorm, xnorm, var
-
- """
- %-----------------------------------------------------------------------
- % End of lsqr.m
- %-----------------------------------------------------------------------
- """
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/dionysus/viewer/__init__.py Tue Jun 05 17:04:10 2012 -0700
@@ -0,0 +1,1 @@
+from diagram import *
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/bindings/python/dionysus/viewer/diagram.py Tue Jun 05 17:04:10 2012 -0700
@@ -0,0 +1,73 @@
+from PyQt4 import QtGui, QtCore
+
+class DiagramPoint(QtGui.QGraphicsEllipseItem):
+ def __init__(self,x,y,radius, p, viewer):
+ super(QtGui.QGraphicsEllipseItem, self).__init__(x,y,radius, radius)
+ self.p = p
+ self.viewer = viewer
+
+ def mousePressEvent(self, event):
+ self.viewer.selection = self.p
+ self.viewer.close()
+
+class DiagramViewer(QtGui.QGraphicsView):
+ def __init__(self, dgm):
+ super(QtGui.QGraphicsView, self).__init__()
+
+ self.selection = None
+
+ self.setRenderHint(QtGui.QPainter.Antialiasing)
+ self.scene = QtGui.QGraphicsScene(self)
+ self.setScene(self.scene)
+
+ minx = min(p[0] for p in dgm)
+ miny = min(p[1] for p in dgm)
+ maxx = max(p[0] for p in dgm)
+ maxy = max(p[1] for p in dgm)
+
+ self.draw_axes(minx,miny,maxx,maxy)
+
+ radius = min(maxx - minx, maxy - miny)/100
+ self.scene.setSceneRect(minx - 10*radius, miny - 10*radius, (maxx - minx) + 20*radius, (maxy - miny) + 20*radius)
+
+ for p in dgm:
+ x,y = p[0],p[1]
+ item = DiagramPoint(x,y,radius, p, self)
+ self.scene.addItem(item)
+
+ # Flip y-axis
+ self.scale(1,-1)
+
+ # Set the correct view
+ rect = self.scene.itemsBoundingRect()
+ self.fitInView(rect, QtCore.Qt.KeepAspectRatio)
+
+ def draw_axes(self, minx, miny, maxx, maxy):
+ # Draw axes and diagonal
+ self.scene.addItem(QtGui.QGraphicsLineItem(0,0, maxx, 0))
+ self.scene.addItem(QtGui.QGraphicsLineItem(minx,0, 0, 0))
+ self.scene.addItem(QtGui.QGraphicsLineItem(0,0, 0, maxy))
+ self.scene.addItem(QtGui.QGraphicsLineItem(0,miny, 0, 0))
+ self.scene.addItem(QtGui.QGraphicsLineItem(0,0, min(maxx, maxy), min(maxx, maxy)))
+ self.scene.addItem(QtGui.QGraphicsLineItem(max(minx,miny), max(minx,miny), 0,0))
+
+ # Dashed, gray integer lattice
+ pen = QtGui.QPen(QtCore.Qt.DashLine)
+ pen.setColor(QtCore.Qt.gray)
+ for i in xrange(int(minx) + 1, int(maxx) + 1):
+ line = QtGui.QGraphicsLineItem(i,0, i, maxy)
+ line.setPen(pen)
+ self.scene.addItem(line)
+ for i in xrange(int(miny) + 1, int(maxy) + 1):
+ line = QtGui.QGraphicsLineItem(0,i, maxx, i)
+ line.setPen(pen)
+ self.scene.addItem(line)
+
+
+def show_diagram(dgm):
+ app = QtGui.QApplication([])
+ view = DiagramViewer(dgm)
+ view.show()
+ view.raise_()
+ app.exec_()
+ return view.selection