pyx/helper.py

   1 #!/usr/bin/env python
   2 # -*- coding: ISO-8859-1 -*-
   3 #
   4 #
   5 # Copyright (C) 2002-2004 Jörg Lehmann <joergl@users.sourceforge.net>
   6 # Copyright (C) 2002-2006 André Wobst <wobsta@users.sourceforge.net>
   7 # Copyright (C) 2005 Michael Schindler <m-schindler@users.sourceforge.net>
   8 #
   9 # This file is part of PyX (http://pyx.sourceforge.net/).
  10 #
  11 # PyX is free software; you can redistribute it and/or modify
  12 # it under the terms of the GNU General Public License as published by
  13 # the Free Software Foundation; either version 2 of the License, or
  14 # (at your option) any later version.
  15 #
  16 # PyX is distributed in the hope that it will be useful,
  17 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19 # GNU General Public License for more details.
  20 #
  21 # You should have received a copy of the GNU General Public License
  22 # along with PyX; if not, write to the Free Software
  23 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  24
  25 import math, types
  26
  27 def sign(x):
  28     return (x >= 0) and 1 or -1
  29
  30 try:
  31     import Numeric, LinearAlgebra
  32     _has_numeric = 1
  33 except:
  34     _has_numeric = 0
  35
  36 def evrealtrisym(d, e):
  37     """returns eigenvalues of a tridiagonal real symmetric matrix
  38
  39     d are the diagonal elements and e are the subdiagonal elements.
  40     len(d) == len(e) + 1; d and e are destroyed
  41     """
  42     n = len(d)
  43     assert len(d) == len(e) + 1
  44     e = e + [0.0]
  45     for l in range(n):
  46         for iter in range(30):
  47             for m in range(l, n-1):
  48                 dd = abs(d[m]) + abs(d[m+1])
  49                 if abs(e[m]) + dd == dd:
  50                     break
  51             else:
  52                 m = n-1
  53             if m == l:
  54                 break
  55             g = (d[l+1]-d[l]) / (2.0*e[l])
  56             r = math.hypot(g, 1.0)
  57             if g >= 0:
  58                 g = d[m]-d[l]+e[l]/(g+abs(r))
  59             else:
  60                 g = d[m]-d[l]+e[l]/(g-abs(r))
  61             s = c = 1.0
  62             p = 0.0
  63             for i in range(m, l, -1):
  64                 f = s*e[i-1]
  65                 b = c*e[i-1]
  66                 e[i] = r = math.hypot(f, g)
  67                 if r == 0.0:
  68                     d[i] -= p
  69                     e[m] = 0.0
  70                     break
  71                 s = f/r
  72                 c = g/r
  73                 g = d[i]-p
  74                 r = (d[i-1]-g)*s + 2.0*c*b
  75                 p = s*r
  76                 d[i] = g+p
  77                 g = c*r-b
  78             else:
  79                 d[l] -= p
  80                 e[l] = g
  81                 e[m] = 0.0
  82         else:
  83             raise RuntimeError("Too many iterations in evrealtrisym")
  84     return d
  85
  86 def realsymtotrisym(a):
  87     """creates a real tridiagonal matrix out of a real symmetric matrix
  88
  89     Eigenvalues are not altered by the transformation; a is the matrix, i.e. a
  90     list of lists but only the diagonal and off-diagonal elements left of the
  91     diagonal are accessed (the other values don't need to be provided and the
  92     lists might be shortend); a is destroyed
  93     """
  94     n = len(a)
  95     d = [None] * n
  96     e = [None] * (n-1)
  97     x = [None] * (n-1)
  98     for i in range(n-1, 0, -1):
  99         l = i
 100         h = scale = 0.0
 101         if l > 1:
 102             for k in range(l):
 103                 scale += abs(a[i][k])
 104             if scale == 0.0:
 105                 e[i-1] = a[i][l-1]
 106             else:
 107                 for k in range(l):
 108                     a[i][k] /= scale
 109                     h += a[i][k]*a[i][k]
 110                 f = a[i][l-1]
 111                 g = math.sqrt(h)
 112                 if f >= 0:
 113                     g = -g
 114                 e[i-1] = scale*g
 115                 h -= f*g
 116                 a[i][l-1] = f-g
 117                 f = 0.0
 118                 for j in range(l):
 119                     g = 0.0
 120                     for k in range(j+1):
 121                         g += a[j][k]*a[i][k]
 122                     for k in range(j+1, l):
 123                         g += a[k][j]*a[i][k]
 124                     x[j] = g/h
 125                     f += x[j]*a[i][j]
 126                 hh = 0.5*f/h
 127                 for j in range(l):
 128                     f = a[i][j]
 129                     x[j] = g = x[j]-hh*f
 130                     for k in range(j+1):
 131                         a[j][k] -= f*x[k] + g*a[i][k]
 132         else:
 133             e[i-1] = a[i][l-1]
 134         d[i] = h
 135     for i in range(n):
 136         d[i] = a[i][i]
 137     return d, e
 138
 139 def evrealsym(a):
 140     """returns eigenvalues of a real symmetric matrix
 141
 142     a is a real symmetric matrix, i.e. a list of lists but only the diagonal
 143     and off-diagonal elements below the diagonal are accessed (other values are
 144     not accessed, i.e. the lists might be shortend).
 145     """
 146     # make a copy and ensure floats
 147     a = [[float(x) for x in row] for row in a]
 148     return evrealtrisym(*realsymtotrisym(a))
 149
 150 def balancerealsquare(a, radix=2):
 151     """modify a general real square matrix to become balanced
 152
 153     This routine replaces the general real square matrix a by a balanced matrix
 154     with identical eigenvalues. A symmetric matrix is already balanced and is
 155     unaffected by this procedure. The parameter radix should be the machine's
 156     floating-point radix to prevent any numerical inaccuracies to be introduced
 157     by this method.
 158     """
 159     n = len(a)
 160     radix2 = radix*radix
 161     last = 0
 162     while not last:
 163         last = 1
 164         for i in range(n):
 165             r = c = 0.0
 166             for j in range(n):
 167                 if j+1 != i+1:
 168                     c += abs(a[j][i])
 169                     r += abs(a[i][j])
 170             if c != 0 and r != 0:
 171                 g = r/radix
 172                 f = 1.0
 173                 s = c+r
 174                 while c < g:
 175                     f *= radix
 176                     c *= radix2
 177                 g = r*radix
 178                 while c > g:
 179                     f /= radix
 180                     c /= radix2
 181                 if (c+r)/f < 0.95*s:
 182                     last = 0
 183                     g = 1.0/f
 184                     for j in range(n):
 185                         a[i][j] *= g
 186                     for j in range(n):
 187                         a[j][i] *= f
 188
 189 def hessrealsquare(a):
 190     """modify a general real square matrix to the upper hessenberg form
 191
 192     This routine replaces the general real square matrix a by a hessenberg
 193     formed matrix with identical eigenvalues.
 194     """
 195     n = len(a)
 196     for m in range(1, n-1):
 197         x = 0.0
 198         i = m
 199         for j in range(m, n):
 200             if (abs(a[j][m-1]) > abs(x)):
 201                 x = a[j][m-1]
 202                 i = j
 203         if i != m:
 204             for j in range(m-1, n):
 205                 a[i][j], a[m][j] = a[m][j], a[i][j]
 206             for j in range(n):
 207                 a[j][i], a[j][m] = a[j][m], a[j][i]
 208         if x:
 209             for i in range(m+1, n):
 210                 y = a[i][m-1]
 211                 if y:
 212                     y /= x
 213                     a[i][m-1] = y
 214                     for j in range(m, n):
 215                         a[i][j] -= y*a[m][j]
 216                     for j in range(n):
 217                         a[j][m] += y*a[j][i]
 218     for i in range(2, n):
 219         for j in range(i-1):
 220             a[i][j] = 0
 221
 222 def evhessrealsquare(a, realonly=0):
 223     """returns eigenvalues of a real square matrix in upper hessenberg form
 224
 225     a is a real square matrix in upper hessenberg form. a is destroyed.
 226     The return value is a list containing a mixture of floats and complex
 227     values.
 228     """
 229     n = len(a)
 230     e = []
 231
 232     anorm = 0.0
 233     for i in range(n):
 234         for j in range(n):
 235             anorm += abs(a[i][j])
 236     t = 0.0
 237     while n > 0:
 238         its = 0
 239         first = 1
 240         while first or l < n-2:
 241             first = 0
 242             for l in range(n-1, 0, -1):
 243                 s = abs(a[l-1][l-1]) + abs(a[l][l])
 244                 if s:
 245                     s = anorm
 246                 if abs(a[l][l-1]) + s == s:
 247                     break
 248             else:
 249                 l = 0
 250             x = a[n-1][n-1]
 251             if l == n-1:
 252                 e.append(x+t)
 253                 n -= 1
 254             else:
 255                 y = a[n-2][n-2]
 256                 w = a[n-1][n-2]*a[n-2][n-1]
 257                 if l == n-2:
 258                     p = 0.5*(y-x)
 259                     q = p*p+w
 260                     z = math.sqrt(abs(q))
 261                     x += t
 262                     if q >= 0.0:
 263                         if p >= 0:
 264                             z = p + abs(z)
 265                         else:
 266                             z = p - abs(z)
 267                         e.append(x+z)
 268                         if z:
 269                             e.append(x-w/z)
 270                         else:
 271                             e.append(x+z)
 272                     elif not realonly:
 273                         e.append(x+p - z*1J)
 274                         e.append(x+p + z*1J)
 275                     n -= 2
 276                 else:
 277                     if its == 30:
 278                         raise RuntimeError("Too many iterations in evhessrealsquare")
 279                     if its == 10 or its == 20:
 280                         t += x
 281                         for i in range(n):
 282                             a[i][i] -= x
 283                         s = abs(a[n-1][n-2]) + abs(a[n-2][n-3])
 284                         y = x = 0.75*s
 285                         w = -0.4375*s*s
 286                     its += 1
 287                     for m in range(n-3, l-1, -1):
 288                         z = a[m][m]
 289                         r = x-z
 290                         s = y-z
 291                         p = (r*s-w)/a[m+1][m] + a[m][m+1]
 292                         q = a[m+1][m+1]-z-r-s
 293                         r = a[m+2][m+1]
 294                         s = abs(p)+abs(q)+abs(r)
 295                         p /= s
 296                         q /= s
 297                         r /= s
 298                         if m == l:
 299                             break
 300                         u = abs(a[m][m-1])*(abs(q)+abs(r))
 301                         v = abs(p) * (abs(a[m-1][m-1]) + abs(z) + abs(a[m+1][m+1]))
 302                         if u + v == v:
 303                             break
 304                     for i in range(m+2, n):
 305                         a[i][i-2] = 0.0
 306                         if i != m+2:
 307                             a[i][i-3]=0.0
 308                     for k in range(m, n-1):
 309                         if k != m:
 310                             p = a[k][k-1]
 311                             q = a[k+1][k-1]
 312                             r = 0.0
 313                             if k != n-2:
 314                                 r = a[k+2][k-1]
 315                             x = abs(p) + abs(q) + abs(r)
 316                             if x:
 317                                 p /= x
 318                                 q /= x
 319                                 r /= x
 320                         if p >= 0:
 321                             s = math.sqrt(p*p+q*q+r*r)
 322                         else:
 323                             s = -math.sqrt(p*p+q*q+r*r)
 324                         if s:
 325                             if k == m:
 326                                 if l != m:
 327                                     a[k][k-1] = -a[k][k-1]
 328                             else:
 329                                 a[k][k-1] = -s*x
 330                             p += s
 331                             x = p/s
 332                             y = q/s
 333                             z = r/s
 334                             q /= p
 335                             r /= p
 336                             for j in range(k, n):
 337                                 p = a[k][j]+q*a[k+1][j]
 338                                 if k != n-2:
 339                                     p += r*a[k+2][j]
 340                                     a[k+2][j] -= p*z
 341                                 a[k+1][j] -= p*y
 342                                 a[k][j] -= p*x
 343                             mmin = min(n-1, k+3)
 344                             for i in range(l, mmin+1):
 345                                 p = x*a[i][k] + y*a[i][k+1]
 346                                 if k != n-2:
 347                                     p += z*a[i][k+2]
 348                                     a[i][k+2] -= p*r
 349                                 a[i][k+1] -= p*q
 350                                 a[i][k] -= p
 351                             else:
 352                                 i = mmin + 1
 353     return e
 354
 355 def evrealsquare(a, copy=1, realonly=0):
 356     """returns eigenvalues of a real square matrix
 357
 358     a is a real square matrix, i.e. a list of lists. The return value is a list
 359     containing a mixture of floats and complex values.
 360     """
 361     if copy:
 362         # make a copy and ensure floats
 363         a = [[float(x) for x in row] for row in a]
 364     balancerealsquare(a)
 365     hessrealsquare(a)
 366     return evhessrealsquare(a, realonly=realonly)
 367
 368 def realpolyroots(coeffs, epsilon=1e-5):
 369     """returns the roots of a polynom with given coefficients
 370
 371     polynomial with coefficients given in coeffs:
 372       0 = \sum_i x^i coeffs[i]
 373     The solution is found via an equivalent eigenvalue problem
 374     """
 375     if not coeffs:
 376         return []
 377     try:
 378         f = -1.0 / coeffs[-1]
 379     except ZeroDivisionError:
 380         return realpolyroots(coeffs[:-1])
 381     else:
 382         n = len(coeffs) - 1
 383         if _has_numeric:
 384             a = Numeric.zeros((n, n), Numeric.Float)
 385             for i in range(n-1):
 386                 a[i+1][i] = 1
 387             for i in range(n):
 388                 a[0][i] = f*coeffs[n-1-i]
 389             roots = []
 390             for root in LinearAlgebra.eigenvalues(a):
 391                 if type(root) == types.ComplexType:
 392                     if not root.imag:
 393                         roots.append(root.real)
 394                 else:
 395                     roots.append(root)
 396             return roots
 397         else:
 398             a = [[0.0 for i in range(n)] for j in range(n)]
 399             for i in range(n-1):
 400                 a[i+1][i] = 1.0
 401             for i in range(n):
 402                 a[0][i] = f*coeffs[n-1-i]
 403             return evrealsquare(a, copy=0, realonly=1)