Changed implementation of WeightedRips to store simplex values (max distance between simplices' vertices) as an invisible layer on top of each simplex object, so that the data() field of WeightedRips has been freed for use by the users again.
# 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
%-----------------------------------------------------------------------
"""