Package sfc :: Package symbolic_utils :: Module symbolic_utils
[hide private]
[frames] | no frames]

Source Code for Module sfc.symbolic_utils.symbolic_utils

  1  #!/usr/bin/env python 
  2  # -*- coding: utf-8 -*- 
  3  """ 
  4  This module contains utility functions for symbolic computations, 
  5  mostly a thin layer on top of swiginac. 
  6  """ 
  7   
  8  # Copyright (C) 2007-2009 Martin Sandve Alnes and Simula Resarch Laboratory 
  9  # 
 10  # This file is part of SyFi. 
 11  # 
 12  # SyFi is free software: you can redistribute it and/or modify 
 13  # it under the terms of the GNU General Public License as published by 
 14  # the Free Software Foundation, either version 2 of the License, or 
 15  # (at your option) any later version. 
 16  # 
 17  # SyFi is distributed in the hope that it will be useful, 
 18  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 19  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 20  # GNU General Public License for more details. 
 21  # 
 22  # You should have received a copy of the GNU General Public License 
 23  # along with SyFi. If not, see <http://www.gnu.org/licenses/>. 
 24  # 
 25  # First added:  2007-10-16 
 26  # Last changed: 2009-03-19 
 27   
 28  import swiginac 
 29  import SyFi 
 30   
 31  from sfc.symbolic_utils.symbol_factory import symbol, symbols, symbolic_vector, symbolic_matrix 
 32   
 33  # ... Utilities for symbolic computations, mostly on matrices 
 34   
 35  _p = symbols( ["x", "y", "z"] ) 
 36   
 37  # ... Some type stuff, shouldn't really be in this file? 
 38   
39 -def is_indexed_type(A):
40 """Checks if the argument is a matrix or lst.""" 41 return isinstance(A, list) \ 42 or isinstance(A, tuple) \ 43 or isinstance(A, swiginac.lst) \ 44 or isinstance(A.evalm(), swiginac.matrix)
45
46 -def is_scalar(e):
47 return not is_indexed_type(e.evalm())
48
49 -def as_matrix(A):
50 """Convert A to a swiginac.matrix from a lst or list.""" 51 if isinstance(A, swiginac.lst): return swiginac.lst_to_matrix(A).evalm() 52 if isinstance(A, list): return swiginac.matrix(len(A), 1, A).evalm() # TODO: handle list of lists? 53 if isinstance(A, tuple): return swiginac.matrix(len(A), 1, A).evalm() # TODO: handle list of lists? 54 if isinstance(A.evalm(), swiginac.matrix): return A.evalm() 55 raise TypeError("ERROR: cannot convert A to a swiginac.matrix, A is of type %s" % str(type(A)))
56
57 -def zeros(m, n):
58 """Returns a m x n matrix with zeros.""" 59 return swiginac.matrix(m, n)
60
61 -def ones(m, n):
62 """Returns a m x n matrix with ones.""" 63 A = swiginac.matrix(m, n) 64 for i in xrange(m): 65 for j in xrange(n): 66 A[i,j] = 1 67 return A
68
69 -def Id(n):
70 """Returns the n x n identity matrix.""" 71 return swiginac.unit_matrix(n)
72
73 -def add(A, B):
74 """Adds two swiginac expressions and calls evalm() on the result before returning.""" 75 if is_indexed_type(A): A = as_matrix(A) 76 if is_indexed_type(B): B = as_matrix(B) 77 return (A+B).evalm()
78
79 -def sub(A, B):
80 """Subtracts two swiginac expressions and calls evalm() on the result before returning.""" 81 if is_indexed_type(A): A = as_matrix(A) 82 if is_indexed_type(B): B = as_matrix(B) 83 return (A-B).evalm()
84
85 -def mul(A, B):
86 """Multiplies two swiginac expressions and calls evalm() on the result before returning.""" 87 if is_indexed_type(A): A = as_matrix(A) 88 if is_indexed_type(B): B = as_matrix(B) 89 return (A*B).evalm()
90
91 -def cross(a, b):
92 """Takes the cross product of two vectors.""" 93 a = as_matrix(a) 94 b = as_matrix(b) 95 if len(a) != 3 or len(b) != 3: 96 raise TypeError("Need 3D vectors a and b in cross(a,b).") 97 im = [a[1], a[2], 98 b[1], b[2]] 99 jm = [a[2], a[0], 100 b[2], b[0]] 101 km = [a[0], a[1], 102 b[0], b[1]] 103 return swiginac.matrix(3, 1, [swiginac.matrix(2,2,m).determinant() for m in (im, jm, km)])
104
105 -def inner(A, B): # FIXME: use multiplication for matrices instead of contraction?
106 """Takes the inner product of A and B: multiplication for scalars, dot product for vectors and lsts, contraction for matrices.""" 107 if is_indexed_type(A): A = as_matrix(A) 108 if is_indexed_type(B): B = as_matrix(B) 109 if isinstance(A, swiginac.matrix) and isinstance(B, swiginac.matrix): 110 if len(A) == len(B): 111 return sum( [A.op(i)*B.op(i) for i in xrange(len(A))] ) 112 if A.cols() == 1: A = A.transpose() # row vector to column vector 113 if B.rows() == 1: B = B.transpose() # column vector to row vector 114 if A.cols() == B.rows(): # matrix-vector or vector-matrix inner product 115 return mul(A, B) 116 raise ValueError("Unmatching size of operands") 117 return (A*B).evalm() 118
119 -def dot(a, b): # FIXME: use multiplication for matrices instead of contraction?
120 return inner(a, b) 121
122 -def contract(A,B):
123 """See inner(A,B).""" 124 return inner(A,B)
125
126 -def det(A):
127 """Returns the determinant of the argument.""" 128 return A.evalm().determinant()
129
130 -def abs(x):
131 """Returns the absolute value of the argument.""" 132 return swiginac.abs(x)
133
134 -def transpose(A):
135 """Returns the transpose of the argument.""" 136 return swiginac.transpose(A.evalm())
137
138 -def trace(A):
139 """Returns the trace of the argument.""" 140 return swiginac.trace(A.evalm())
141
142 -def inverse(A):
143 """Returns the inverse of the argument.""" 144 A = A.evalm() 145 if isinstance(A, swiginac.matrix): 146 return A.evalm().inverse() 147 return 1.0 / A
148
149 -def rank(x):
150 if isinstance(x, swiginac.matrix): 151 r = x.rows() 152 c = x.cols() 153 if r > 1 and c > 1: return 2 154 if r*c > 1: return 1 155 return 0
156
157 -def shape(x):
158 if isinstance(x, swiginac.matrix): 159 r = x.rows() 160 c = x.cols() 161 if r > 1 and c > 1: return (r, c) 162 if r*c > 1: return (r*c,) 163 return (1,)
164
165 -def diff(f, x):
166 """Returns df/dx, where x is a symbol or symbolic matrix.""" 167 if type(x) == swiginac.matrix: 168 f_rank = rank(f) 169 x_rank = rank(x) 170 if f_rank+x_rank > 2: raise RuntimeError("Cannot apply diff to f,x with ranks %d,%d" % (f_rank, x_rank)) 171 # FIXME: support different combinations of ranks for f,x 172 A = zeros(x.rows(), x.cols()) 173 for i in range(A.rows()): 174 for j in range(A.cols()): 175 sym = x[i,j] 176 assert isinstance(sym, swiginac.symbol) 177 val = f # [...] FIXME 178 A[i,j] = swiginac.diff(val, sym) 179 return A 180 assert isinstance(x, swiginac.symbol) 181 return swiginac.diff(f, x)
182
183 -def ddx(f, i, GinvT=None):
184 """Returns df/dx_i, where i is the number of the coordinate. 185 GinvT is an optional geometry mapping.""" 186 # without mapping: 187 if GinvT is None: 188 return diff(f, x(i)) 189 # with mapping: 190 dfdx = 0 191 for k in range(GinvT.rows()): 192 dfdx += GinvT[i, k] * swiginac.diff(f, _p[k]) 193 return dfdx
194 195
196 -def grad(u, GinvT=None):
197 """Returns the gradient of u w.r.t. the spacial variables x,y[,z], depending on the dimension specified by SyFi.initSyFi(nsd).""" 198 199 ran = rank(u) 200 sh = shape(u) 201 202 if ran == 0: 203 n = SyFi.cvar.nsd 204 r, c = n, 1 205 u = swiginac.matrix(1,1,[u]) 206 elif ran == 1: 207 n = sh[0] 208 r, c = n, n 209 else: 210 raise ValueError("Taking gradient of a matrix is not supported.") 211 212 if GinvT is None: 213 GinvT = Id(n) 214 215 g = zeros(r, c) 216 for i in range(r): 217 for j in range(c): 218 f = u[j,0] 219 for k in range(n): 220 g[i,j] += GinvT[i,k] * diff(f, _p[k]) 221 return g
222 223
224 -def div(u, GinvT=None):
225 r = u.rows() 226 c = u.cols() 227 ran = rank(u) 228 if ran == 2: 229 if r != c: 230 raise RuntimeError("Taking divergence of an unsymmetric matrix.") 231 n = r 232 if GinvT is None: 233 GinvT = Id(n) 234 d = zeros(r, 1) 235 for j in range(r): 236 for i in range(r): 237 for k in range(r): 238 d[j,0] += GinvT[i,k] * diff(u[i,j], _p[k]) 239 elif ran == 1: 240 n = r*c 241 if GinvT is None: 242 GinvT = Id(n) 243 d = 0 244 for i in range(n): 245 for k in range(n): 246 d += GinvT[i,k] * diff(u[i,0], _p[k]) 247 else: 248 raise RuntimeError("Taking divergence of a scalar.") 249 return d
250 251
252 -def curl(u, GinvT=None):
253 g = grad(u, GinvT) 254 c = swiginac.matrix(2,1) 255 c[0,0] = -g[1,0] 256 c[1,0] = g[0,0] 257 return c
258 259
260 -def laplace(u, GinvT=None):
261 return div(grad(u, GinvT), GinvT)
262 263 264 265 if __name__ == '__main__': 266 a = swiginac.matrix(3,1,[1,0,0]) 267 b = swiginac.matrix(3,1,[0,1,0]) 268 c = swiginac.matrix(3,1,[0,0,1]) 269 270 print "" 271 print cross(a,b) 272 print cross(a,c) 273 print cross(c,b) 274 print cross(b,c) 275 276 import SyFi 277 SyFi.initSyFi(2) 278 279 x, y, z = symbols(["x", "y", "z"]) 280 x2, y2, z2 = symbols(["x", "y", "z"]) 281 assert (x-x2).is_zero() 282 assert (y-y2).is_zero() 283 assert (z-z2).is_zero() 284 285 p2 = swiginac.matrix(2, 1, [x,y]) 286 p3 = swiginac.matrix(3, 1, [x,y,z]) 287 288 289 A = swiginac.matrix(2,2, [x*x*x, x*x*y, x*y*y, y*y*y]) 290 sh = shape(A) 291 assert rank(A) == 2 292 assert len(sh) == 2 293 assert sh[0] == 2 294 assert sh[1] == 2 295 print "" 296 print "A", A 297 print "div A", div(A) 298 299 300 v = swiginac.matrix(2, 1, [3*x,5*y]) 301 sh = shape(v) 302 assert rank(v) == 1 303 assert len(sh) == 1 304 assert sh[0] == 2 305 print "" 306 print "v", v 307 print "div v", div(v) 308 print "grad v", grad(v) 309 print "dv/dx ", diff(v, p2) 310 311 312 v = swiginac.matrix(2, 1, [3*y,5*x]) 313 sh = shape(v) 314 assert rank(v) == 1 315 assert len(sh) == 1 316 assert sh[0] == 2 317 print "" 318 print "v", v 319 print "div v", div(v) 320 print "grad v", grad(v) 321 print "dv/dx ", diff(v, p2) 322