1from __future__ import division 2import numpy as np 3import numpy.matlib as ml 4import numpy.linalg as la 5import scipy.linalg as sla 6import scipy.optimize as opt 7import scipy.signal as sig 8import sympy as sp 9import ctypes as ct 10import slicot 11import mpmath 12import collections, warnings, tempfile, os, operator, sys, shutil, time 13import commands, subprocess, math, re, logging 14 15import dk_templates, generate_code 16from models import GND, Out, Node 17from dk_lib import printoptions 18 19def pkgconfig(*packages, **kw): 20 flag_map = {'-I': 'include_dirs', '-L': 'library_dirs', '-l': 'libraries'} 21 for token in commands.getoutput("pkg-config --libs --cflags %s" % ' '.join(packages)).split(): 22 kw.setdefault(flag_map.get(token[:2]), []).append(token[2:]) 23 return kw 24 25def matrix_add(*args): 26 l = [m for m in args if m.size] 27 if l: 28 return reduce(ml.add, l) 29 else: 30 return args[0] 31 32 33class TransformOptions(object): 34 35 def __init__(self, **kw): 36 self._set(**kw) 37 38 def _set(self, partition=False, chain=False, svd_prec=None, decompose=False): 39 self.partition = partition 40 self.chain = chain 41 self.svd_prec = svd_prec 42 self.decompose = decompose 43 44 def copy(self, **kw): 45 t = TransformOptions() 46 t.__dict__.update(self.__dict__) 47 t._set(kw) 48 return t 49 50 @property 51 def svd_prec(self): 52 if self._svd_prec is None: 53 self._svd_prec = math.sqrt(sys.float_info.epsilon) 54 return self._svd_prec 55 56 @svd_prec.setter 57 def svd_prec(self, v): 58 self._svd_prec = v 59 60 61class ConvergenceError(Exception): pass 62 63class CodeWrapError(Exception): pass 64 65class CodeWrapper(object): 66 67 _module_counter = 0 68 69 def __init__(self, d, code, tempdir=None, verbose=False, flags=[]): 70 self.code = code 71 if tempdir is not None: 72 tempdir = os.path.abspath(tempdir) 73 self.filepath = tempdir 74 self.flags = flags 75 self.quiet = not verbose 76 packages = ["eigen3"] 77 if d["method"] in ("hybr", "lm"): 78 packages.append('cminpack') 79 dct = pkgconfig(*packages) 80 libs = dct.get("libraries",[]) 81 libs.append("m") 82 if d["resample"]: 83 libs.append("zita-resampler") 84 self.script_dict = dict( 85 includes = " ".join("-I%s" % v for v in dct.get("include_dirs",[])), 86 libraries = " ".join("-l%s" % v for v in set(libs)), 87 #defines = "-DCHECK_BOUNDS", 88 defines = "", 89 debug = False, 90 sourcename = "%s_%d.cpp" % (d["id"], CodeWrapper._module_counter), 91 soname = "%s_%d.so" % (d["id"], CodeWrapper._module_counter), 92 soname_ = "%s_%d_.so" % (d["id"], CodeWrapper._module_counter), 93 optimize = False,#(d["method"] == "table"), 94 ) 95 self.build_script = "./build_script_%d" % CodeWrapper._module_counter 96 97 def _process_files(self, routine): 98 command = self.command 99 command.extend(self.flags) 100 null = open(os.devnull, 'w') 101 try: 102 if self.quiet: 103 retcode = subprocess.call(command, stdout=null, stderr=subprocess.STDOUT) 104 else: 105 retcode = subprocess.call(command) 106 except OSError: 107 retcode = 1 108 if retcode: 109 raise CodeWrapError( 110 "Error while executing command: %s" % " ".join(command)) 111 112 @property 113 def command(self): 114 command = [self.build_script] 115 return command 116 117 def _prepare_files(self): 118 with open(self.build_script, 'w') as f: 119 print >> f, dk_templates.build_script_template.render(self.script_dict) 120 os.fchmod(f.fileno(), 0o777) 121 122 123 def _generate_code(self): 124 for fname, content in self.code.items(): 125 if fname == "c_source": 126 fname = self.script_dict['sourcename'] 127 with open(fname, 'w') as f: 128 f.write(content) 129 130 def wrap_code(self, load_module=True): 131 workdir = self.filepath or tempfile.mkdtemp("_dk_compile") 132 if not os.access(workdir, os.F_OK): 133 os.mkdir(workdir) 134 oldwork = os.getcwd() 135 os.chdir(workdir) 136 try: 137 self._generate_code() 138 self._prepare_files() 139 self._process_files(None) 140 path = os.path.join(workdir, self.script_dict["soname"]) 141 if load_module: 142 return path, SimulateC(path) 143 else: 144 return path 145 finally: 146 CodeWrapper._module_counter +=1 147 os.chdir(oldwork) 148 if not self.filepath: 149 shutil.rmtree(workdir) 150 151 152class LinearFilter(object): 153 154 def __init__(self, p, jacobi): 155 self.in_mat = p.N["I"] 156 self.out_mat = p.N["O"] 157 if self.in_mat.shape[0] != 1 or self.out_mat.shape[0] != 1: 158 raise ValueError("linear filter generation only implemented for circuits with one input and one output channel") 159 pot_func = p.get_pot_funcs() 160 Rv = sp.diag(*[1/(f * v) for (a, f), v in zip(pot_func, p.Pv)]) 161 Nv = p.N["P"] 162 S = p.S + Nv.T * Rv * Nv 163 if jacobi is not None: 164 S -= p.N["Nr"].T * jacobi * p.N["Nl"] 165 self.S = S 166 ss = set() 167 self.syms = [] 168 for a, f in pot_func: 169 if a not in ss: 170 self.syms.append(a) 171 ss.add(a) 172 self.subst_var_default = p.get_variable_defaults() 173 174 def get_s_coeffs(self): 175 s = sp.symbols("s") 176 expr = self.solve(self.S, self.in_mat, self.out_mat) 177 r = [sp.poly(e, s) for e in sp.fraction(expr)] 178 tc = r[1].TC() 179 for v in tc.atoms(sp.Symbol): 180 tc = tc.subs(v, 1) 181 b_coeffs, b_terms = self.collect_s_coeffs(sp.poly(r[0], s), 'b', tc) 182 a_coeffs, a_terms = self.collect_s_coeffs(sp.poly(r[1], s), 'a', tc) 183 if a_coeffs[0] == -1: 184 a_coeffs *= -1 185 b_coeffs *= -1 186 return b_coeffs, a_coeffs, b_terms / a_terms 187 188 def collect_s_coeffs(self, expr, prefix, tc): 189 s = sp.symbols("s") 190 monoms = expr.monoms() 191 max_degree = monoms[0][0] 192 filter_coeffs = np.zeros(max_degree+1, dtype=object) 193 ll = 0 194 for e, i in zip(expr.coeffs(), monoms): 195 i = i[0] 196 if self.syms: 197 # factorize according to variable symbols 198 x = sp.poly(e, self.syms) 199 ss = 0 200 for co, o in zip(x.coeffs(), x.monoms()): 201 ss += reduce(operator.mul, [pow(y, p) for y, p in zip(self.syms, o)], 1) * (co/tc).simplify() 202 else: 203 ss = e/tc 204 filter_coeffs[i] = ss 205 sym = "%s%d" % (prefix, i) 206 ll += sp.symbols(sym) * pow(s, i) 207 return filter_coeffs, ll 208 209 def transform_bilinear(self, expr): 210 s = sp.symbols("s") 211 fs = sp.symbols("fs") 212 c = sp.symbols("c") 213 z = sp.symbols("z") 214 b = c*(z-1)/(z+1) 215 r = [sp.poly(e, z) for e in sp.fraction(expr.subs(s, b).ratsimp())] 216 return self.collect_z_symbolic(r[0]), self.collect_z_symbolic(r[1]), 2 * fs 217 218 def collect_z_symbolic(self, expr): 219 c = sp.symbols("c") 220 monoms = expr.monoms() 221 max_degree = monoms[0][0] 222 filter_coeffs = np.zeros(max_degree+1, dtype=object) 223 for e, i in zip(expr.coeffs(), monoms): 224 idx = max_degree - i[0] 225 filter_coeffs[idx] = e 226 return filter_coeffs 227 228 def convert_variable_dict(self, subst_var): 229 d = self.subst_var_default.copy() 230 df = set(subst_var.keys()) - set([str(v) for v in d.keys()]) 231 if df: 232 raise ValueError("unknown variable(s): %s" % ", ".join(df)) 233 d.update(dict([(sp.symbols(k), v) for k, v in subst_var.items()])) 234 return d 235 236 def get_z_coeffs(self, samplerate=None, subst_var=None, as_expr=True): 237 if samplerate is None: 238 c = 2 * sp.symbols("fs") 239 else: 240 c = 2 * samplerate 241 z = sp.symbols("z") 242 b = c*(z-1)/(z+1) 243 expr = self.S.subs(sp.symbols("s"), b) 244 expr = self.solve(expr, self.in_mat, self.out_mat) 245 if subst_var is not None: 246 expr = expr.subs(subst_var) 247 syms = None 248 else: 249 syms = self.syms 250 # divide coeffs by magnitude of trailing denominator coefficient == a1 coefficient of filter 251 lc = sp.poly(sp.fraction(expr)[1], z).LC() 252 lc = lc.subs(self.subst_var_default) 253 if samplerate is None: 254 lc = lc.subs(sp.symbols("fs"), 48000.) 255 r = [sp.poly((e/lc).expand(), z) for e in sp.fraction(expr)] 256 return (self.collect_z_coeffs(r[0], 'b', syms, as_expr), 257 self.collect_z_coeffs(r[1], 'a', syms, as_expr)) 258 259 def collect_z_coeffs(self, expr, prefix, syms, as_expr=True): 260 monoms = expr.monoms() 261 if as_expr: 262 max_degree = monoms[0][0] 263 filter_coeffs = np.zeros(max_degree+1, dtype=object) 264 else: 265 l = [] 266 for e, i in zip(expr.coeffs(), monoms): 267 i = i[0] 268 if syms: 269 # factorize according to variable symbols 270 x = sp.poly(e.expand(), syms) 271 if not as_expr: 272 l.append(zip(x.monoms(), x.coeffs())) 273 else: 274 ss = 0 275 for co, o in zip(x.coeffs(), x.monoms()): 276 ss += reduce(operator.mul, [pow(y, p) for y, p in zip(syms, o)], 1) * co.evalf() 277 else: 278 if not as_expr: 279 l.append(((0,), e)) 280 else: 281 ss = e.evalf() 282 if as_expr: 283 idx = max_degree-i 284 filter_coeffs[idx] = sp.horner(ss, wrt=(syms or [])+[sp.symbols("fs")]) 285 if as_expr: 286 return filter_coeffs 287 else: 288 return l, syms 289 290 def coeffs_as_faust_code(self, prefix, coeffs): 291 l = [] 292 for i, c in enumerate(coeffs): 293 c = re.sub(r'([a-zA-Z0-9]+)\*\*(\d+)', r'pow(\1,\2)', str(c)) 294 l.append('%s%d = %s;' % (prefix, i, c)) 295 return l 296 297 def print_coeffs(self, prefix, coeffs, f=sys.stdout): 298 print >>f, "\n\n".join(self.coeffs_as_faust_code(prefix, coeffs)) 299 300 def spectrum(self, b, a, start_freq=20, stop_freq=10000, fs=48000, nbins=8*1024): 301 b = [float(x) for x in b] 302 a = [float(x) for x in a] 303 w, h = sig.freqz(b, a, nbins) 304 cut = slice(int(round(2 * nbins * start_freq / fs)), int(round(2 * nbins * stop_freq / fs)) + 1) 305 w = (w[cut] * (fs/(2*np.pi))) 306 h = 20*np.log10(abs(h[cut])) 307 return w, h 308 309 def solve(self, S, in_mat, out_mat): 310 v = sp.Matrix(sp.symbols("v:%d" % self.S.shape[0])) 311 try: 312 p = subprocess.Popen(("maxima","-b","/dev/fd/0","--very-quiet","2>&1 >/dev/null"),shell=True, 313 stdin=subprocess.PIPE, stdout=subprocess.PIPE, 314 stderr=subprocess.PIPE) 315 except OSError as e: 316 raise RuntimeError( 317 "can't start maxima -- please check maxima installation [%s]" % e) 318 out, expr = p.communicate( 319 "stringout(\"/dev/fd/2\",facsum(linsolve([%s], [%s])[%d], s));\n" % ( 320 ", ".join(["%s = %s" % e for e in zip(S * v, in_mat)]), 321 ", ".join([str(sym) for sym in v]), 322 np.array(out_mat).nonzero()[1]+1, 323 )) 324 def maxima_error(s): 325 print out 326 print "----" 327 print "Truncated output: " + expr[:100] + "..." 328 raise ValueError(s) 329 syms = set() 330 for i in S: 331 syms |= i.atoms(sp.Symbol) 332 e = expr.split("=",1) 333 if len(e) != 2: 334 maxima_error("unexpected output from maxima") 335 e = e[1].rstrip(";\n").replace("^","**") 336 try: 337 return eval(e, dict([(str(sym),sym) for sym in syms])) 338 except Exception as ex: 339 maxima_error("can't eval maxima output [%s]" % ex) 340 341def get_state_transform_trace(A, B, C): 342 "return the trace of the gram matrices (sensitivity measure)" 343 Wc = ml.matrix(sla.solve_discrete_lyapunov(A.A, (B*B.T).A)) 344 Wo = ml.matrix(sla.solve_discrete_lyapunov(A.T.A, (C.T*C).A)) 345 #R = ml.matrix(sla.cholesky(Wo, lower=False)) 346 #X = R * Wc * R.T 347 #U, s, V = sla.svd(X) 348 #U = ml.matrix(U) 349 #V = ml.matrix(V) 350 #lli = np.sqrt(np.sqrt(s)) 351 #U = R.I * V.T * ml.diag(lli) 352 #Ui = ml.diag(1/lli) * V * R 353 #return U, Ui 354 return np.trace(Wc), np.trace(Wo) 355 356def get_state_transform(A, B, C, tol=0): 357 "calculate balanced reduced state space model realization" 358 # balance matrices 359 A, B, C, maxred, scale = slicot.tb01id('A', A, B, C) 360 # transform to schur form 361 A, B, C, U, WR, WI = slicot.tb01wd(A, B, C) 362 # reduce and balance 363 T, Ti, A, B, C, hsv = slicot.ab09ax('D', 'B', A, B, C, tol=tol) 364 # return transformation matrix (+ inverse) and hankel singular values 365 return (np.matrix(np.diag(scale).dot(U).dot(T)), 366 np.matrix(Ti.dot(U.T).dot(np.diag(1/scale))), 367 hsv) 368 369 370class NonlinEquations(object): 371 372 def __init__(self, eq, v_slice): 373 self.eq = eq 374 self.v_slice = v_slice 375 self.p_slice = v_slice 376 self.i_slice = v_slice 377 self.nn = self.v_slice.stop - self.v_slice.start 378 self.nni = self.nn 379 self.nno = self.nn 380 self.U = self.Mi = ml.matrix(np.identity(self.nn)) 381 self.Hc = ml.zeros((self.nn, 1)) 382 self.pins = eq.nlin_elements[v_slice] 383 self.name = "%s:%d" % (Parser.format_element(np.sort(self.pins, axis=0)[0]), self.nn) 384 self.subblocks = [] 385 386 @staticmethod 387 def create(eq, K, CZ, v_slice, opts): 388 "return a Permution and an instance of NonlinEquations or derived" 389 Pn = np.arange(len(CZ)) 390 # permute nonlinear part to make left upper submatrix of K blockdiagonal 391 p, blocklist = NonlinEquations.get_block_indices(K[v_slice][:,v_slice], CZ[v_slice]) 392 if len(blocklist) > 1: 393 Pn[v_slice] = p 394 K = K[Pn][:,Pn] 395 CZ = CZ[Pn] 396 if opts.partition and len(blocklist) > 1: 397 nlin = PartitionedNonlinEquations(eq, v_slice) 398 Pn = Pn[nlin.create(eq, K, CZ, blocklist, opts)] 399 return Pn, nlin 400 # permute the system if there is more than one strongly connected component 401 p, blocklist = NonlinEquations.find_scc(K[v_slice][:,v_slice], eq.get_parser()) 402 if opts.chain and len(blocklist) > 1: 403 Pn[v_slice] = p 404 K = K[Pn][:,Pn] 405 CZ = CZ[Pn] 406 nlin = ChainedNonlinEquations(eq, v_slice) 407 Pn = Pn[nlin.create(eq, K, CZ, blocklist, opts)] 408 return Pn, nlin 409 return np.arange(len(CZ)), NonlinEquations(eq, v_slice) 410 411 @staticmethod 412 def rebase_nonlinear_functions(f, Pn): 413 Pni = np.argsort(Pn) # inverse permutation 414 a = np.zeros(len(f), dtype=object) 415 for j, (expr, vl, base) in enumerate(f[Pn]): 416 a[j] = (expr, vl, Pni[base]) 417 return a 418 419 @staticmethod 420 def get_unique_rows(a): 421 """returns unique rows of a 2-dim array 422 """ 423 if not a.size: 424 return a 425 order = np.lexsort(a.T) 426 a = a[order] 427 diff = np.diff(a, axis=0) 428 ui = np.ones(len(a), dtype=bool) 429 ui[1:] = (diff != 0).any(axis=1) 430 return a[ui] 431 432 @staticmethod 433 def get_blockdiag_permutation(B): 434 """returns the permuted index list to transform B to block diagonal form 435 436 B: quadratic (sparse) matrix 437 returns: index permutation (p), block indexlist (bl) 438 439 To get the block diagonal matrix A: 440 A = B[p][:,p] 441 442 Permutation matrix: 443 Q = matrix(eye(len(p)))[p] 444 A = Q * B * Q.T 445 446 To permute a list of row or column labels with the matrix: 447 labels_for_A = [labels_for_B[i] for i in p] 448 449 To get the ith block matrix: 450 r = slice(bl[i], bl[i+1]) 451 A[r,r] 452 """ 453 B = ml.matrix(B, bool) 454 G = B.T * B 455 while True: 456 G1 = G * G 457 if (G1 == G).all(): 458 break 459 G = G1 460 n = B.shape[0] 461 p = [] 462 bl = [0] 463 for row in NonlinEquations.get_unique_rows(G.A): 464 i = np.nonzero(row)[0] 465 if len(i): 466 p.extend(i) 467 bl.append(len(p)) 468 if len(p) != B.shape[0]: 469 p.extend(set(range(B.shape[0])) - set(p)) 470 return p, bl 471 472 @staticmethod 473 def get_block_indices(K, CZ): 474 "returns a permutation and a list of component slices" 475 nz = len(CZ) 476 n = np.count_nonzero(CZ) 477 if n != nz: 478 pc = np.argsort(CZ == 0) 479 K = K[pc][:,pc][:n,:n] 480 p, bl = NonlinEquations.get_blockdiag_permutation(K) 481 if n != nz: 482 p = np.array(pc)[p+range(n, nz)] 483 return p, [slice(i,j) for i, j in zip(bl[:-1], bl[1:])] 484 485 @staticmethod 486 def decompose(a): 487 U, s, V = la.svd(a, full_matrices=False) 488 o = np.sqrt(np.sum(s*s) / np.count_nonzero(s)) 489 sr = np.where(s, s, o) 490 V1 = ml.diag(sr) * V 491 V1i = (ml.diag(1/sr) * V).T 492 U1 = U * np.diag(np.where(s, 1, 0)) 493 if (U1 <= 0).all(): 494 U1 = -U1 495 V1 = -V1 496 V1i = -V1i 497 return U1, V1, V1i 498 499 @staticmethod 500 def find_scc(K, parser): 501 from scipy.sparse import csr_matrix 502 from scipy.sparse.csgraph import connected_components 503 n, label = connected_components(csr_matrix(K), connection="strong") 504 V = np.zeros(n) 505 L = [[] for i in range(n)] 506 for i, l in enumerate(L): 507 col = K[:,label==i] 508 for j in range(n): 509 if i != j and col[label==j].any(): 510 l.append(j) 511 V[j] += 1 512 p = [] 513 l = [] 514 N = np.arange(len(K)) 515 while True: 516 idx = np.argwhere(V == 0) 517 if not idx.size: 518 break 519 V[idx] -= 1 520 for i in idx: 521 pi = N[label==i] 522 j = len(p) 523 l.append(slice(j, j+len(pi))) 524 p.extend(pi) 525 V[L[int(i)]] -= 1 526 return p, l 527 528 def make_input_trans(self, matrix_list, opts): 529 if self.eq.np or opts.svd_prec < 0: 530 return matrix_list 531 # find number of linear independent inputs to the nonlinear function 532 M = np.concatenate(matrix_list, axis=1) 533 if not M.size: 534 return matrix_list 535 U, SV, V = la.svd(M, full_matrices=False) 536 if SV[0] != 0: 537 SV = np.where(SV / SV[0] > opts.svd_prec, SV, 0) 538 nni = np.count_nonzero(SV) 539 if nni == self.nn: 540 return matrix_list 541 U = U[:,:nni] 542 UH = U.H 543 self.U = U 544 self.nni = nni 545 self.p_slice = slice(self.p_slice.start, self.p_slice.start+nni) 546 return [UH * m for m in matrix_list] 547 548 def make_output_trans(self, matrix_list, opts): 549 if self.nn == 0 or self.eq.np != 0 or opts.svd_prec < 0: 550 return matrix_list 551 # find number of linear independent outputs of the nonlinear function 552 M = np.concatenate(matrix_list, axis=0) 553 U, SV, V = la.svd(M, full_matrices=False) 554 if SV[0] != 0: 555 SV = np.where(SV / SV[0] > opts.svd_prec, SV, 0) 556 n = np.count_nonzero(SV) 557 if n < self.nn and n < self.eq.nx + self.eq.no: 558 self.nno = n 559 self.Mi = (np.diag(SV) * V)[:n] 560 s = np.cumsum([0]+[m.shape[0] for m in matrix_list]) 561 matrix_list = [U[slice(*sl),:n] for sl in zip(s[:-1],s[1:])] 562 elif self.eq.nx + self.eq.no < self.nn: 563 self.nno = self.eq.nx + self.eq.no 564 self.Mi = M 565 M = np.eye(self.nno, self.nno) 566 s = np.cumsum([0]+[m.shape[0] for m in matrix_list]) 567 matrix_list = [M[slice(*sl)] for sl in zip(s[:-1],s[1:])] 568 self.i_slice = slice(self.i_slice.start, self.i_slice.start+self.nno) 569 return matrix_list 570 571 def transform_p0(self, p0, K, v0): 572 return p0 573 574 575class ChainedNonlinEquations(NonlinEquations): 576 577 def create(self, eq, K, CZ, blocklist, opts): 578 self.subblocks = [] 579 Pn = np.arange(len(K)) 580 for sl in blocklist: 581 p, nlin = NonlinEquations.create(eq, K, CZ, sl, opts.copy(partition=False)) 582 self.subblocks.append(nlin) 583 Pn = Pn[p] 584 return Pn 585 586 def make_input_trans(self, matrix_list, opts): 587 if opts.svd_prec < 0: 588 return matrix_list 589 i = self.subblocks[0].p_slice.start 590 j = self.subblocks[-1].p_slice.stop 591 ll = [[] for m in matrix_list] 592 if i: 593 for l, m in zip(ll, matrix_list): 594 l.append(m[0:i]) 595 for bl in self.subblocks: 596 r = bl.make_input_trans([m[bl.p_slice] for m in matrix_list], opts) 597 for l, m in zip(ll, r): 598 l.append(m) 599 for l, m in zip(ll, matrix_list): 600 if j < m.shape[0]: 601 l.append(m[j:]) 602 return [np.concatenate(l, axis=0) for l in ll] 603 604 def make_output_trans(self, matrix_list, opts): 605 if opts.svd_prec < 0: 606 return matrix_list 607 F, C = matrix_list ##FIXME 608 i = self.subblocks[0].i_slice.start 609 j = self.subblocks[-1].i_slice.stop 610 l_F = [] 611 l_C = [] 612 if i: 613 l_F.append(F[:,0:i]) 614 l_C.append(C[:,0:i]) 615 for bl in self.subblocks: 616 Fs, Cs = bl.make_output_trans((F[:,bl.i_slice], C[:,bl.i_slice]), opts) 617 l_F.append(Fs) 618 l_C.append(Cs) 619 if j < F.shape[0]: 620 l_F.append(F[:,j:]) 621 l_C.append(C[:,j:]) 622 return np.concatenate(l_F, axis=1), np.concatenate(l_C, axis=1) 623 624 625class PartitionedNonlinEquations(NonlinEquations): 626 627 def create(self, eq, K, CZ, blocklist, opts): 628 self.cc_slice = slice(blocklist[-1].stop, self.v_slice.stop) 629 self.subblocks = [] 630 Pn = np.arange(len(K)) 631 for sl in blocklist: 632 p, nlin = NonlinEquations.create(eq, K, CZ, sl, opts.copy(chain=False)) 633 self.subblocks.append(nlin) 634 Pn = Pn[p] 635 return Pn 636 637 def transform_p0(self, p0, K, v0): 638 endv = self.cc_slice.start 639 endp = self.p_slice.stop - (self.cc_slice.stop - self.cc_slice.start) 640 p = p0.copy() 641 p[:endp] = (p0 + self.Hc)[:endp] + self.Ku * np.matrix(v0[endv:]).T 642 return p 643 644 def make_input_trans(self, matrix_list, opts): 645 if opts.svd_prec < 0: 646 Ku = [self.eq.K[bl.p_slice,self.cc_slice] for bl in self.subblocks] 647 self.Ku = np.concatenate(Ku, axis=0) 648 return matrix_list 649 l = [] 650 j = 0 651 Ku = [] 652 for bl in self.subblocks: 653 mlist = [m[bl.p_slice] for m in matrix_list] 654 mlist.append(self.eq.K[bl.p_slice,self.cc_slice]) 655 bl.Hc = self.Hc[bl.p_slice] 656 mlist = bl.make_input_trans(mlist, opts) 657 bl.p_slice = slice(j, bl.p_slice.stop - (bl.p_slice.start - j)) 658 j = bl.p_slice.stop 659 l.append(mlist[:-1]) 660 Ku.append(mlist[-1]) 661 self.Ku = np.concatenate(Ku, axis=0) 662 Hc = self.Hc[self.cc_slice] 663 l.append([m[self.cc_slice] for m in matrix_list]) 664 n_old = self.p_slice.stop - self.p_slice.start 665 cc = self.cc_slice.stop - self.cc_slice.start 666 self.p_slice = slice(self.p_slice.start, j+cc) 667 self.nni = self.p_slice.stop - self.p_slice.start 668 self.Hc = ml.zeros((self.nni, 1)) 669 self.Hc[j:] = Hc 670 self.U = ml.eye(n_old, self.nni, self.nni-n_old) 671 j = 0 672 for bl in self.subblocks: 673 self.U[j:j+bl.U.shape[0],bl.p_slice] = bl.U 674 j += bl.U.shape[0] 675 return [np.concatenate(m, axis=0) for m in zip(*l)] 676 677 def make_output_trans(self, matrix_list, opts): 678 if opts.svd_prec < 0: 679 Kl = [self.eq.K[self.cc_slice, bl.i_slice] for bl in self.subblocks] 680 Kl.append(self.eq.K[self.cc_slice, self.cc_slice]) 681 self.Kl = np.concatenate(Kl, axis=1) 682 return matrix_list 683 l = [] 684 j = 0 685 Kl = [] 686 for bl in self.subblocks: 687 mlist = [m[:,bl.i_slice] for m in matrix_list] 688 mlist.append(self.eq.K[self.cc_slice, bl.i_slice]) 689 mlist = bl.make_output_trans(mlist, opts) 690 bl.i_slice = slice(j, bl.i_slice.stop - (bl.i_slice.start - j)) 691 j = bl.i_slice.stop 692 l.append(mlist[:-1]) 693 Kl.append(mlist[-1]) 694 Kl.append(self.eq.K[self.cc_slice, self.cc_slice]) 695 self.Kl = np.concatenate(Kl, axis=1) 696 l.append([m[:,self.cc_slice] for m in matrix_list]) 697 n_old = self.i_slice.stop - self.i_slice.start 698 cc = self.cc_slice.stop - self.cc_slice.start 699 self.i_slice = slice(self.i_slice.start, j+cc) 700 self.nno = self.i_slice.stop - self.i_slice.start 701 self.Mi = ml.eye(self.nno, n_old, n_old-self.nno) 702 j = 0 703 for bl in self.subblocks: 704 self.Mi[bl.i_slice, j:j+bl.Mi.shape[1]] = bl.Mi 705 j += bl.Mi.shape[1] 706 return [np.concatenate(m, axis=1) for m in zip(*l)] 707 708 def decompose_blocks(self, F, C): 709 Kni = ml.identity(self.nn) 710 end = self.subblocks[-1].v_slice.stop 711 #self.Kl = self.eq.K[end:, :].copy() 712 #return F, C 713 for bl in self.subblocks: 714 sl = bl.v_slice 715 U, V, Vi = self.decompose(self.eq.K[end:, sl]) 716 m = ml.zeros_like(self.Kl[:,sl]) 717 m[:,:U.shape[1]] = U 718 self.Kl[:,sl] = m 719 #self.Mi[sl,sl] = V 720 bl.Mi = V 721 m = ml.zeros_like(Kni[sl,sl]) 722 m[:,:Vi.shape[1]] = Vi 723 Kni[sl,sl] = m 724 return F * Kni, C * Kni 725 726 727class EquationSystem(object): 728 729 def __init__(self, parser, jacobi_par=None, opts=None): 730 self.parser = parser 731 self.jacobi_par = jacobi_par 732 if opts is None: 733 opts = TransformOptions() 734 735 self.nx, self.nn, self.ni, self.no, self.np = [len(parser.element_name[j]) for j in 'X', 'N', 'I', 'O', 'P'] 736 self.nlin_elements = np.array(parser.element_name["N"]) 737 self.Nxl, self.Nxr, self.Nnl, self.Nnr, self.No, self.Nv, self.I = [ 738 parser.N[j] for j in "Xl", "Xr", "Nl", "Nr", "O", "P", "I"] 739 self.CV = parser.ConstVoltages 740 m = parser.mm 741 self.f = parser.get_nonlin_funcs() 742 self.CZ = parser.CZ 743 if self.jacobi_par is not None: 744 J, Jc, select, unselect = self.jacobi_par 745 self.S = parser.S - self.Nnr[unselect].T * J[unselect][:,unselect] * self.Nnl[unselect] 746 self.CV = self.CV + (self.Nnr[unselect].T * Jc[unselect]).T 747 self.Nnr = self.Nnr[select] 748 self.Nnl = self.Nnl[select] 749 self.nn = len(select) 750 self.CZ = self.CZ[select] 751 self.f = NonlinEquations.rebase_nonlinear_functions(self.f, select) 752 parser.f = self.f 753 parser.element_name["N"] = np.array(parser.element_name["N"])[select] 754 else: 755 self.S = parser.S 756 self.Si = self.S.I 757 758 if self.nn: 759 K = (self.Nnl * self.Si * self.Nnr.T != 0) 760 for j, (expr, vl, base) in enumerate(self.f): 761 for i in base: 762 K[j, i] = True 763 v_slice = slice(0, len(K)) 764 Pn, self.nonlin = NonlinEquations.create(self, K, self.CZ, v_slice, opts) 765 self.apply_permutation(Pn) 766 else: 767 self.nonlin = None 768 769 Z = ml.diag(parser.Z) 770 self.Tx = m * Z * self.Nxl * self.Si 771 self.A = self.Tx * self.Nxr.T - (Z if parser.TR else ml.diag((parser.Z - 1) / 2.0)) 772 self.B = self.Tx * self.I.T 773 self.Bc = self.Tx * self.CV.T 774 self.C = self.Tx * self.Nnr.T 775 776 self.To = self.No * self.Si 777 self.D = self.To * self.Nxr.T 778 self.E = self.To * self.I.T 779 self.Ec = self.To * self.CV.T 780 self.F = self.To * self.Nnr.T 781 782 self.Tn = self.Nnl * self.Si 783 self.G = self.Tn * self.Nxr.T 784 self.H = self.Tn * self.I.T 785 self.Hc = self.Tn * self.CV.T 786 self.K = self.Tn * self.Nnr.T 787 788 if self.nn: 789 self.nonlin.Hc = self.Hc 790 self.G, self.H = self.nonlin.make_input_trans((self.G, self.H), opts) 791 self.F, self.C = self.nonlin.make_output_trans((self.F, self.C), opts) 792 793 if self.np: 794 # Woodbury Identity: \left(A+UCV \right)^{-1} = A^{-1} - A^{-1}U \left(C^{-1}+VA^{-1}U \right)^{-1} VA^{-1}, 795 self.Twl = self.Si * self.Nv.T 796 self.Q = self.Nv * self.Twl 797 self.Uxl = m * Z * self.Nxl * self.Twl 798 self.Uo = self.No * self.Twl 799 self.Unl = self.Nnl * self.Twl 800 self.Twr = self.Si.T * self.Nv.T 801 self.Uxr = self.Nxr * self.Twr 802 self.Uu = self.I * self.Twr 803 self.Ucv = self.CV * self.Twr 804 self.Unr = self.Nnr * self.Twr 805 806 self.mp_cols = self.nx + self.ni 807 808 def apply_permutation(self, Pn): 809 self.nlin_elements = self.nlin_elements[Pn] 810 self.CZ = self.CZ[Pn] 811 self.f = NonlinEquations.rebase_nonlinear_functions(self.f, Pn) 812 self.Nnl = self.Nnl[Pn] 813 self.Nnr = self.Nnr[Pn] 814 815 def transform_state(self, Ts, Tsi): 816 self.A = Tsi * self.A * Ts 817 self.B = Tsi * self.B 818 self.Bc = Tsi * self.Bc 819 self.C = Tsi * self.C 820 self.D = self.D * Ts 821 self.G = self.G * Ts 822 self.nx = self.A.shape[0] 823 self.mp_cols = self.nx + self.ni 824 if self.np: 825 self.Uxl = Tsi * self.Uxl 826 self.Uxr = Ts.T * self.Uxr 827 828 def get_parser(self): 829 return self.parser 830 831 def get_Mo(self): 832 return np.concatenate((self.D, self.E, self.F), axis=1) 833 834 def get_Mx(self): 835 return np.concatenate((self.A, self.B, self.C), axis=1) 836 837 def get_Mp(self): 838 return np.concatenate((self.G, self.H), axis=1) 839 840 def get_UR(self): 841 assert hasattr(self, "Q") 842 return np.concatenate((self.Uxr.T, self.Uu.T, self.Unr.T), axis=1) 843 844 def get_mx_cols(self): 845 return self.A.shape[1] + self.B.shape[1] + self.C.shape[1] 846 847 def get_mp_cols(self): 848 return self.G.shape[1] + self.H.shape[1] 849 850 def get_npl(self): 851 return len(self.parser.pot_list) 852 853 def get_status(self): 854 return "nx=%d, nni=%d, nn=%d, nno=%d, ni=%d, no=%d, np=%d" % (self.nx, self.nni, self.nn, self.nno, self.ni, self.no, self.np) 855 856 857class Simulate(object): 858 859 def __init__(self, eq, solver): 860 self.eq = eq 861 self.parser = eq.get_parser() 862 if solver is None: 863 self.solver_dict = dict() 864 elif isinstance(solver, basestring): 865 self.solver_dict = dict(method = solver) 866 else: 867 self.solver_dict = dict(solver) 868 self.solver_method = self.solver_dict.get("method", 'hybr') 869 self.max_homotopy_iter = self.solver_dict.get("max_homotopy_iter", 1000) 870 try: 871 del self.solver_dict["max_homotopy_iter"] 872 except KeyError: 873 pass 874 try: 875 del self.solver_dict["method"] 876 except KeyError: 877 pass 878 879 def get_solver(self): 880 d = self.solver_dict.copy() 881 d["method"] = self.solver_method 882 d["max_homotopy_iter"] = self.max_homotopy_iter 883 return d 884 885 def get_eq(self): 886 return self.eq 887 888 def get_parser(self): 889 return self.parser 890 891 892class SimulatePy(Simulate): 893 894 def __init__(self, eq, solver=None, dc_method="A"): 895 Simulate.__init__(self, eq, solver) 896 self.dc_method = dc_method 897 parser = eq.get_parser() 898 self.pot_func = parser.get_pot_funcs() 899 self.Pv = parser.Pv 900 self.pot = parser.pot 901 self.out_labels = parser.out_labels() 902 self.v0 = parser.V.get("v0", np.zeros(eq.nn)) 903 self.x = np.zeros(eq.nx) 904 self.o0 = np.zeros(eq.no) 905 self.op = parser.op 906 self.compile_py_func() 907 self.calc_dc(parser.op, method=dc_method) 908 909 def set_variable(self, var, val): 910 self.pot[var] = val 911 912 def calc_Rv(self): 913 n = self.eq.np 914 Rv = ml.matrix(np.zeros((n, n))) 915 for i, ((a, f), p) in enumerate(zip(self.pot_func, self.Pv)): 916 k = str(a) 917 try: 918 v = self.pot[k] 919 except KeyError: 920 v = self.pot[k] = 0.5 921 Rv[i, i] = float(f.subs({a: v})) * p 922 return Rv 923 924 def calc_matrixes(self): 925 eq = self.eq 926 if eq.np: 927 Qi = (eq.Q + self.calc_Rv()).I 928 Tx = eq.Uxl * Qi 929 To = eq.Uo * Qi 930 Tn = eq.Unl * Qi 931 return (eq.A - Tx * eq.Uxr.T, 932 eq.B - Tx * eq.Uu.T, 933 eq.Bc - Tx * eq.Ucv.T, 934 eq.C - Tx * eq.Unr.T, 935 eq.D - To * eq.Uxr.T, 936 eq.E - To * eq.Uu.T, 937 eq.Ec - To * eq.Ucv.T, 938 eq.F - To * eq.Unr.T, 939 eq.G - Tn * eq.Uxr.T, 940 eq.H - Tn * eq.Uu.T, 941 eq.Hc - Tn * eq.Ucv.T, 942 eq.K - Tn * eq.Unr.T, 943 ) 944 else: 945 return (eq.A, eq.B, eq.Bc, eq.C, 946 eq.D, eq.E, eq.Ec, eq.F, 947 eq.G, eq.H, eq.Hc, eq.K, 948 ) 949 950 def calc_linear_sys_matrices(self, v0=None): 951 eq = self.eq 952 if eq.nn: 953 nonlin = eq.nonlin 954 J = self.jacobi(v0) 955 X = nonlin.Mi * (np.diag(eq.CZ) - J * eq.K).I * J * nonlin.U 956 # equivalent formula for X: 957 # X = nonlin.Mi * J * (np.diag(eq.CZ) - eq.K * J).I * nonlin.U 958 A = eq.A + eq.C * X * eq.G 959 B = eq.B + eq.C * X * eq.H 960 D = eq.D + eq.F * X * eq.G 961 return A, B, D 962 else: 963 return eq.A, eq.B, eq.D 964 965 def balance_realization(self, tol=None): 966 if self.eq.np: 967 return 968 A, B, D = self.calc_linear_sys_matrices() 969 T, Ti, hsv = get_state_transform(A, B, D, tol=tol) 970 logger = logging.getLogger("balance") 971 if T.shape[0] != T.shape[1]: 972 logger.info("reduced system size %d -> %d" % T.shape) 973 logger.debug("HSV = %s" % hsv) 974 self.eq.transform_state(T, Ti) 975 self.calc_dc(self.eq.parser.op, method=self.dc_method) 976 977 def solve(self, func, v0, args=(), method='hybr', options=None): 978 if method == "lm" and options and "maxfev" in options: 979 options["maxiter"] = options["maxfev"] 980 del options["maxfev"] 981 try: 982 with warnings.catch_warnings(): 983 warnings.filterwarnings(action="error") 984 res = opt.root(func, v0, args=args, method=method, options=options) 985 except RuntimeWarning as e: 986 raise ConvergenceError(e) 987 else: 988 if not res.success: 989 raise ConvergenceError(res.message) 990 return res 991 992 def solve_using_homotopy(self, func, v0, method='hybr', options=None): 993 # use homotopy 994 points = [0, 1] 995 max_iter = 100 996 for tries in range(max_iter): 997 try: 998 res = self.solve(func, v0, args=(points[1],), method=method, options=options) 999 except ConvergenceError as e: 1000 msg = e 1001 points.insert(1, (points[0]+points[1])/2) 1002 continue 1003 if len(points) == 2: 1004 return res 1005 v0 = res.x 1006 points = points[1:] 1007 raise ConvergenceError("more than %d iterations (list msg: %s)" % (max_iter, msg)) 1008 1009 def calc_dc(self, u, method="A"): 1010 u = ml.matrix(u, dtype=np.float64).T 1011 A, B, Bc, C, D, E, Ec, F, G, H, Hc, K = self.calc_matrixes() 1012 if A.size == 0: 1013 if len(self.eq.f): 1014 p = self.eq.nonlin.U * H * u 1015 def func(v, fact): 1016 return (p + Hc * fact + K * self.calc_i(v) - ml.matrix(self.eq.CZ * v).T).A1 1017 self.v0 = self.solve_using_homotopy(func, self.v0).x 1018 else: 1019 I = ml.eye(len(A)) 1020 if method == "A": 1021 try: 1022 Ai = (I - A).I 1023 except la.LinAlgError: 1024 method = "N" 1025 if method == "N": 1026 if self.eq.nonlin is None: 1027 Bu = B * u + Bc 1028 def func(v, fact): 1029 vv = ml.matrix(v).T 1030 return ((A - I) * vv + Bu*fact).A1 1031 res = self.solve_using_homotopy(func, self.x) 1032 self.x = res.x 1033 else: 1034 G1 = ml.append(self.eq.nonlin.U * G, A, axis=0) 1035 H1 = ml.append(self.eq.nonlin.U * H, B, axis=0) 1036 K1 = ml.append(K, C * self.eq.nonlin.Mi, axis=0) 1037 Hc1 = ml.append(Hc, Bc, axis=0) 1038 CZ1 = np.append(self.eq.CZ, np.ones(len(self.x))) 1039 n = len(self.v0) 1040 def func(v, fact): 1041 return (G1 * ml.matrix(v[n:]).T + H1 * u + Hc1*fact + K1 * self.calc_i(v) - ml.matrix(CZ1 * v).T).A1 1042 res = self.solve_using_homotopy(func, np.append(self.v0, self.x)) 1043 self.v0 = res.x[:n] 1044 self.x = res.x[n:] 1045 else: 1046 self.x = matrix_add(Ai * B * u, Ai * Bc) 1047 if K.size != 0: 1048 T = self.eq.nonlin.U * G * Ai 1049 p = matrix_add(T * B * u, self.eq.nonlin.U * H * u) 1050 Hc1 = matrix_add(T * Bc, Hc) 1051 KK = T * C * self.eq.nonlin.Mi + K 1052 def func(v, fact): 1053 return (p + Hc1 * fact + KK * self.calc_i(v) - ml.matrix(self.eq.CZ * v).T).A1 1054 self.v0 = self.solve_using_homotopy(func, self.v0).x 1055 self.x += Ai * C * self.eq.nonlin.Mi * self.calc_i(self.v0) 1056 self.x = self.x.A1 1057 self.v00 = self.v0.copy() 1058 self.x0 = self.x.copy() 1059 self.p0 = matrix_add(G * np.matrix(self.x0).T, H * np.matrix(self.op).T) 1060 if self.eq.nonlin: 1061 self.p0 = self.eq.nonlin.transform_p0(self.p0, K, self.v0) 1062 self.o0 = matrix_add(D * np.matrix(self.x0).T, E * np.matrix(self.op).T, Ec, F * self.eq.nonlin.Mi * self.calc_i(self.v0)).A1 1063 self.last_p = self.eq.nonlin.U * self.p0 + Hc 1064 else: 1065 self.o0 = matrix_add(D * np.matrix(self.x0).T, E * np.matrix(self.op).T, Ec).A1 1066 1067 def calc_i(self, v): 1068 i = ml.zeros(len(self.ff)).T 1069 for n, (f, base) in enumerate(self.ff): 1070 i[n] = f(*v[base]) 1071 return i 1072 1073 def solve_one(self, p, K, s, sm=None, sd=None): 1074 if sm is None: 1075 sm = self.solver_method 1076 if sd is None: 1077 sd = self.solver_dict 1078 i = ml.zeros((p.shape[0],1)) 1079 vv = self.v0.copy() 1080 def func(v): 1081 for n, (f, base) in enumerate(self.ff[s]): 1082 ##FIXME 1083 vv[s] = v 1084 i[n] = f(*vv[base]) 1085 return (p + K * i - ml.matrix(v).T).A1 1086 self.v0[s] = self.solve(func, self.v0[s], method=sm, options=sd).x 1087 return i 1088 1089 def nonlin_py(self, p, K, Hc): 1090 p = self.eq.nonlin.U * p + Hc 1091 i = ml.zeros((self.eq.nn, 1)) 1092 if isinstance(self.eq.nonlin, PartitionedNonlinEquations): 1093 rc = slice(self.eq.blocklist[-1].stop, None) 1094 #sm = ["hybr", "lm", "hybr"] 1095 #sd = [dict(), dict(diag=(1e3,1),factor=1e-1), dict()] 1096 #sm = ["lm", "lm", "lm"] 1097 #sd = [dict(diag=(1e3,1),factor=1e-1), dict(diag=(1e3,1),factor=1e-1), dict(diag=(1e3,1),factor=1e-1)] 1098 sm = [None,None,None] 1099 sd = [None,None,None] 1100 def func(icc, fact): 1101 #print "*", fact, icc 1102 p1 = self.last_p + (p - self.last_p) * fact 1103 for j, s in enumerate(self.eq.blocklist): 1104 k = K[s] 1105 try: 1106 i[s] = self.eq.Kn[j] * self.solve_one(p1[s]+k[:,rc].dot(icc).T, k[:,s], s, sm[j], sd[j]) 1107 except ConvergenceError as e: 1108 print "conv error in %d: %s" % (j, e) 1109 raise 1110 #raise SystemExit 1111 i[rc] = icc.reshape(len(icc), 1) 1112 return (p1[rc] + self.eq.Kl * i).A1 1113 self.v0[rc] = self.solve_using_homotopy(func, self.v0[rc], method=self.solver_method, options=self.solver_dict).x 1114 elif isinstance(self.eq.nonlin, ChainedNonlinEquations): 1115 for j, s in enumerate(self.eq.nonlin.blocklist): 1116 pp = p[s] 1117 if s.start: 1118 pp += K[s][:,:s.start] * i[:s.start] 1119 i[s] = self.solve_one(pp, K[s][:,s], s) 1120 i = self.eq.nonlin.Mi * i 1121 else: 1122 def func(v, fact): 1123 i[:] = self.calc_i(v) 1124 p1 = self.last_p + (p - self.last_p) * fact 1125 return (p1 + K * i - ml.matrix(self.eq.CZ * v).T).A1 1126 self.v0 = self.solve_using_homotopy(func, self.v0, method=self.solver_method, options=self.solver_dict).x 1127 i = self.eq.nonlin.Mi * i 1128 self.last_p = p 1129 return i 1130 1131 def compile_py_func(self): 1132 self.ff = np.zeros(len(self.eq.f), dtype=object) 1133 if not self.eq.nn: 1134 # model is linearized 1135 return 1136 for j, (expr, vl, base) in enumerate(self.eq.f): 1137 self.ff[j] = (sp.lambdify(vl, expr, modules=["mpmath", "math", "sympy"]), base) 1138 1139 def calc_di(self, v, j): 1140 i = np.zeros(len(self.eq.f)) 1141 for n, (f, vl, base) in enumerate(self.eq.f): 1142 for k, var in zip(base, vl): 1143 if k == j: 1144 s = dict([(sym,val) for sym, val in zip(vl, v[base])]) 1145 i[n] = f.diff(var).subs(s) 1146 break 1147 else: 1148 i[n] = 0. 1149 return i 1150 1151 def jacobi_numeric(self, v0=None): 1152 if v0 is None: 1153 v0 = self.v00 1154 J = np.zeros((len(self.eq.f), len(v0))) 1155 i0 = self.calc_i(v0).A1 1156 eps = 1e-4 1157 for j in range(len(v0)): 1158 v = v0.copy() 1159 v[j] += eps 1160 J[:, j] = (self.calc_i(v).A1 - i0) / eps 1161 return ml.matrix(J) 1162 1163 def jacobi_symbolic(self, v0=None): 1164 if v0 is None: 1165 v0 = self.v00 1166 J = np.zeros((len(self.eq.f), len(v0))) 1167 for j in range(len(v0)): 1168 J[:, j] = self.calc_di(v0, j) 1169 return ml.matrix(J) 1170 1171 jacobi = jacobi_symbolic 1172 #jacobi = jacobi_numeric 1173 1174 def eval_py(self, v_in, ii=-1): 1175 self.x = ml.matrix(self.x).T ##FIXME 1176 v_in = ml.matrix(v_in) 1177 assert v_in.shape[1] == self.eq.ni 1178 y = np.empty((v_in.shape[0], self.eq.no)) 1179 t1 = time.time() 1180 A, B, Bc, C, D, E, Ec, F, G, H, Hc, K = self.calc_matrixes() 1181 self.minmax = np.array(((float("inf"), float("-inf")),) * G.shape[0]) 1182 for n, u in enumerate(v_in): 1183 u = ml.matrix(u).T 1184 if len(self.v0) == 0: 1185 i = ml.matrix(()).T 1186 else: 1187 try: 1188 p = G * self.x + H * u 1189 self.minmax[:,0] = np.min((self.minmax[:,0], np.ravel(p)), axis=0) 1190 self.minmax[:,1] = np.max((self.minmax[:,1], np.ravel(p)), axis=0) 1191 i = self.nonlin_py(p, K, Hc) 1192 except ConvergenceError as e: 1193 print "##", n 1194 raise 1195 y[n,:] = matrix_add(D * self.x, E * u, Ec, F * i).A1 1196 self.x = matrix_add(A * self.x, B * u, Bc, C * i) 1197 self.time_per_sample = (time.time() - t1)/(n+1) 1198 self.x = self.x.A1 ##FIXME 1199 return y 1200 1201 def __call__(self, v_in, ii=-1): 1202 return self.eval_py(v_in, ii) 1203 1204 1205class CheckedDict(dict): 1206 def __setitem__(self, n, v): 1207 assert n not in self 1208 return dict.__setitem__(self, n, v) 1209 def overwrite(self, n, v): 1210 return dict.__setitem__(self, n, v) 1211 1212class PluginDef(object): 1213 1214 def __init__(self, id): 1215 self.id = id 1216 self.name = id 1217 self._description = None 1218 self.category = "External" 1219 self._shortname = None 1220 self.namespace = id 1221 self.lv2_plugin_type = None 1222 self.lv2_versioned_id = id 1223 self.lv2_minor_version = 0 1224 self.lv2_micro_version = 0 1225 1226 @staticmethod 1227 def _lfmt(s): 1228 if s is None: 1229 return 0 1230 if not s: 1231 return '""' 1232 return 'N_("%s")' % s 1233 1234 @property 1235 def description(self): 1236 return self._description or "" 1237 @description.setter 1238 def description(self, v): 1239 self._description = v 1240 1241 @property 1242 def shortname(self): 1243 return self._shortname or "" 1244 @shortname.setter 1245 def shortname(self, v): 1246 self._shortname = v 1247 1248 @property 1249 def s_id(self): 1250 return '"%s"' % self.id 1251 1252 @property 1253 def l_name(self): 1254 return self._lfmt(self.name) 1255 1256 @property 1257 def l_description(self): 1258 return self._lfmt(self._description) 1259 1260 @property 1261 def l_category(self): 1262 return self._lfmt(self.category) 1263 1264 @property 1265 def l_shortname(self): 1266 return self._lfmt(self._shortname) 1267 1268 1269class BuildCModule(Simulate): 1270 1271 def __init__(self, name, sim, solver=None, c_tempdir=None, c_verbose=False, 1272 c_real="double", extra_sources=None, linearize=False, pre_filter=None, 1273 post_filter=None, generator=None): 1274 if solver is None: 1275 solver = sim.get_solver() 1276 Simulate.__init__(self, sim.get_eq(), solver) 1277 self.name = name 1278 self.c_tempdir = c_tempdir 1279 self.c_verbose = c_verbose 1280 self.c_real = c_real 1281 self.extra_sources = extra_sources 1282 self.pre_filter = pre_filter 1283 self.post_filter = post_filter 1284 if generator is None: 1285 generator = generate_code.CodeGenerator 1286 self.generator = generator 1287 parser = sim.get_parser() 1288 if linearize: 1289 l = [i for i, e in enumerate(parser.element_name["N"]) if "linearize" not in parser.V[e[0]]] 1290 li = [i for i, e in enumerate(parser.element_name["N"]) if "linearize" in parser.V[e[0]]] 1291 J = sim.jacobi() 1292 Jc = sim.calc_i(sim.v00) 1293 par = J, Jc, l, li 1294 self.eq = EquationSystem(parser, par) 1295 sim = SimulatePy(self.eq, solver) 1296 self.v0, self.x0, self.p0, self.o0 = sim.v0, sim.x0, sim.p0, sim.o0 1297 self.op = parser.op 1298 self.pot_func = parser.get_pot_funcs() 1299 self.pot = parser.pot 1300 self.Pv = parser.Pv 1301 self.pot_attr = parser.get_pot_attr() 1302 self.out_labels = parser.out_labels() 1303 self.fs = parser.fs 1304 self.max_homotopy_iter = 64000 1305 self.resample = False 1306 self.solver_params = None 1307 self.dev_interface = True 1308 self.build_script = None 1309 self.plugindef = PluginDef(name) 1310 1311 def get_executor(self): 1312 return self.compile_c_func().wrap_code() 1313 1314 def compile_c_func(self): 1315 eq = self.eq 1316 d = CheckedDict(name=self.name, comment=time.ctime()) 1317 d["solver_maptype"] = "unsigned short" ##FIXME 1318 d["fs"] = self.fs 1319 d["resample"] = self.resample 1320 d["c_real"] = self.c_real 1321 for j in "nx", "nn", "ni", "no", "np", "mp_cols": 1322 d[j] = getattr(eq, j) 1323 for j in "nni", "nno": 1324 d[j] = getattr(eq.nonlin, j, 0) 1325 d["npl"] = self.eq.get_npl() 1326 d["v0_data"] = ",".join([str(j) for j in self.v0]) 1327 d["x0_data"] = ",".join([str(j) for j in self.x0]) 1328 d["p0_data"] = ",".join([str(j) for j in self.p0.A1]) 1329 d["o0_data"] = ",".join([str(j) for j in self.o0]) 1330 d["op_data"] = ",".join([str(j) for j in self.op]) 1331 d["out_labels"] = ",".join(['"%s"' % j for j in self.out_labels]) 1332 d["nlin_elements"] = [Parser.format_element(v) for v in self.eq.nlin_elements] 1333 d["method"] = method = "linear" if eq.nn == 0 else self.solver_method 1334 pot_list = self.eq.get_parser().pot_list 1335 d["pot_vars"] = ",".join(['"%s"' % v for v in pot_list]) 1336 d["pot"] = ",".join([str(self.pot.get(v,0.5)) for v in pot_list]) 1337 d["pre_filter"] = self.pre_filter or "" 1338 d["post_filter"] = self.post_filter or "" 1339 d["post_process"] = "" 1340 d['id'] = d["name"] 1341 d['plugindef'] = self.plugindef 1342 d['build_script'] = self.build_script 1343 d["dev_interface"] = self.dev_interface 1344 solver = self.solver_dict.copy() 1345 solver["method"] = method 1346 solver["max_homotopy_iter"] = self.max_homotopy_iter 1347 if method == "table": 1348 d["extra_sources"] = self.extra_sources 1349 else: 1350 d["extra_sources"] = "" 1351 cg = self.generator( 1352 self.eq, solver, self.solver_params, self.pot, pot_list, 1353 self.pot_func, self.pot_attr, self.Pv, self.extra_sources 1354 ).generate(d) 1355 return CodeWrapper(d, cg, self.c_tempdir, self.c_verbose) 1356 1357 1358class SimulateC(object): 1359 1360 def __init__(self, soname): 1361 self.soname = soname 1362 self.load_from_shared_lib() 1363 self.v0, self.x, self.p0, self.o0, self.op = self.c_get_dc() 1364 self.v00 = self.v0 1365 self.x0 = self.x 1366 self.c_set_state(self.v0, self.x) 1367 self.c_calc_pot_update(np.array([self.pot[v] for v in self.pot_list], dtype=self.dtp)) 1368 self.eval = self.eval_c 1369 self.nonlin = self.nonlin_c 1370 1371 def set_variable(self, var, val): 1372 assert var in self.pot 1373 self.pot[var] = val 1374 self.c_calc_pot_update(np.array([self.pot[v] for v in self.pot_list], dtype=self.dtp)) 1375 1376 def nonlin_c(self, p):#, K, Hc): ##FIXME 1377 return self.c_calc_nonlin(p, self.v0) ##FIXME 1378 1379 def eval_c(self, v_in, ii=-1): 1380 v_in = np.array(v_in) 1381 if v_in.ndim == 1: 1382 v_in = v_in.reshape((len(v_in),1)) 1383 assert v_in.shape[1] == self.ni 1384 y = np.empty((v_in.shape[0], self.no)) 1385 x = self.x 1386 v = self.v0 1387 t1 = time.time() 1388 for n, u in enumerate(v_in): 1389 y[n,:], x, v = self.c_calc(u, x, v) 1390 self.time_per_sample = (time.time() - t1)/n 1391 self.x = x 1392 self.v0 = v 1393 return y 1394 1395 def reset(self): 1396 self.v0 = self.v00 1397 self.x = self.x0 1398 self.c_set_state(self.v00, self.x0) 1399 1400 def eval_c(self, v_in, ii=-1): 1401 self.c_set_state(self.v0, self.x) 1402 t1 = time.time() 1403 try: 1404 return self.c_calc_stream(v_in, ii) 1405 finally: 1406 self.time_per_sample = (time.time() - t1)/v_in.shape[0] 1407 self.v0, self.x, self.minmax, info, nfev, fnorm = self.c_get_info() 1408 1409 def load_from_shared_lib(self): 1410 INTERFACE_VERSION = 5 1411 try: 1412 lib = ct.cdll.LoadLibrary(self.soname) 1413 except OSError as e: 1414 raise SystemExit(e) 1415 try: 1416 version = lib.get_interface_version() 1417 except AttributeError: 1418 raise SystemExit("%s: bad shared lib, missing get_interface_version" % self.soname) 1419 if version != INTERFACE_VERSION: 1420 raise SystemExit("interface version %d expected (found: %d)" % (INTERFACE_VERSION, version)) 1421 c_char_pp = ct.POINTER(ct.c_char_p) 1422 c_get_structure = lib.get_structure 1423 c_get_structure.restype = None 1424 c_get_structure.argtypes = [c_char_pp, ct.POINTER(ct.c_int), ct.POINTER(ct.c_int), ct.POINTER(ct.POINTER(ct.c_int)), 1425 ct.POINTER(ct.POINTER(ct.c_int)), ct.POINTER(c_char_pp), ct.POINTER(c_char_pp), 1426 c_char_pp, ct.POINTER(c_char_pp), ct.POINTER(ct.POINTER(ct.c_double)), 1427 ct.POINTER(c_char_pp), c_char_pp] 1428 nm = ct.c_char_p() 1429 sz = ct.c_int() 1430 fs = ct.c_int() 1431 t = ct.POINTER(ct.c_int)() 1432 tc = ct.POINTER(ct.c_int)() 1433 pins = c_char_pp() 1434 cnm = c_char_pp() 1435 p = ct.c_char_p() 1436 plist = c_char_pp() 1437 pvals = ct.POINTER(ct.c_double)() 1438 ol = c_char_pp() 1439 cmt = ct.c_char_p() 1440 c_get_structure(ct.byref(nm), ct.byref(sz), ct.byref(fs), ct.byref(t), ct.byref(tc), ct.byref(pins), 1441 ct.byref(cnm), ct.byref(p), ct.byref(plist), ct.byref(pvals), ct.byref(ol), 1442 ct.byref(cmt)) 1443 nx, ni, no, npl, nn, nni, nno, nc, end = [t[i] for i in range(9)] 1444 if end != -1: 1445 raise SystemExit("%s: bad sequence length in get_structure") 1446 self.nx, self.ni, self.no, self.npl, self.nn, self.nni, self.nno = nx, ni, no, npl, nn, nni, nno 1447 self.comp_sz = comp_sz = np.array([[slice(tc[6*i+2*j],tc[6*i+2*j+1]) for j in range(3)] for i in range(nc)]) 1448 self.pins = np.array([pins[i] for i in range(nn)]) 1449 self.comp_names = [cnm[2*i] for i in range(nc)] 1450 self.comp_namespace = [cnm[2*i+1] for i in range(nc)] 1451 self.name = nm.value 1452 self.data_size = sz.value 1453 self.fs = fs.value 1454 self.method = p.value 1455 self.pot_list = [plist[i] for i in range(npl)] 1456 self.pot = dict([(plist[i], pvals[i]) for i in range(npl)]) 1457 self.out_labels = [ol[i] for i in range(no)] 1458 self.comment = cmt.value 1459 if self.data_size == 8: 1460 c_real = ct.c_double 1461 self.dtp = np.float64 1462 elif self.data_size == 4: 1463 c_real = ct.c_float 1464 self.dtp = np.float32 1465 else: 1466 raise ValueError("unknown data size: %d" % self.data_size) 1467 def c_arr(n, w=False): 1468 flags = ['C'] 1469 if w: 1470 flags.append('W') 1471 return np.ctypeslib.ndpointer(dtype=c_real, ndim=1, shape=(n,), flags=flags) 1472 def c_mat(w=False): 1473 flags = ['C'] 1474 if w: 1475 flags.append('W') 1476 return np.ctypeslib.ndpointer(dtype=c_real, ndim=2, flags=flags) 1477 c_int_p = ct.POINTER(ct.c_int) 1478 c_real_p = ct.POINTER(c_real) 1479 c_calc = lib.calc 1480 c_calc.restype = ct.c_int 1481 c_calc.argtypes = [c_arr(ni), c_arr(nx), c_arr(nn,True), c_arr(nx,True), c_arr(no,True), c_int_p, c_int_p, c_real_p] 1482 c_set_state = lib.set_state 1483 c_set_state.restype = None 1484 c_set_state.argtypes = [c_arr(nn), c_arr(nx)] 1485 c_get_info = lib.get_info 1486 c_get_info.restype = None 1487 c_get_info.argtypes = [c_arr(nn,True), c_arr(nx, True), c_arr(nni,True), c_arr(nni,True), c_int_p, c_int_p, c_real_p] 1488 c_calc_stream = lib.calc_stream 1489 c_calc_stream.restype = ct.c_int 1490 c_calc_stream.argtypes = [c_mat(), c_mat(True), ct.c_int] 1491 info = ct.c_int() 1492 nfev = ct.c_int() 1493 fnorm = c_real() 1494 if nn: 1495 c_calc_nonlin = lib.calc_nonlin 1496 c_calc_nonlin.restype = ct.c_int 1497 c_calc_nonlin.argtypes = [ct.c_int, c_mat(), c_mat(True), c_arr(nn, True), c_int_p, c_int_p, c_real_p] 1498 def calc_nonlin(p, v): 1499 assert p.shape[1] == nni+npl, (p.shape[1], nni+npl) 1500 p = np.array(p, dtype=self.dtp, order='C', copy=False) 1501 i = np.zeros((p.shape[0], nno), dtype=self.dtp, order='C') 1502 r = c_calc_nonlin(len(p), p, i, v, ct.byref(info), ct.byref(nfev), ct.byref(fnorm)) 1503 if r != 0: 1504 raise ValueError("convergence error: info=%d, nfev=%d, fnorm=%g" % (info.value, nfev.value, fnorm.value)) 1505 return i 1506 self.c_calc_nonlin = calc_nonlin 1507 def define_func(f, i, npl): 1508 v_slice, p_slice, i_slice = comp_sz[i] 1509 nni = p_slice.stop - p_slice.start 1510 nno = i_slice.stop - i_slice.start 1511 dtp = self.dtp 1512 def func(p, v=None): 1513 assert p.shape[1] == nni+npl, (p.shape[1], nni+npl) 1514 if v is None: 1515 v = self.v0 1516 p = np.array(p, dtype=dtp, order='C', copy=False) 1517 i = np.zeros((p.shape[0], nno), dtype=dtp, order='C') 1518 r = f(len(p), p, i, v, ct.byref(info), ct.byref(nfev), ct.byref(fnorm)) 1519 if r != 0: 1520 raise ValueError("convergence error: info=%d, nfev=%d, fnorm=%g" % (info.value, nfev.value, fnorm.value)) 1521 return i 1522 func.name = self.comp_names[i] 1523 func.pins = self.pins[v_slice] 1524 func.v_slice = v_slice 1525 func.p_slice = p_slice 1526 func.i_slice = i_slice 1527 return func 1528 self.c_calc_comp = [] 1529 for i in range(nc): 1530 t = getattr(lib, "calc_%s" % self.comp_namespace[i]) 1531 t.restype = ct.c_int 1532 t.argtypes = [ct.c_int, c_mat(), c_mat(True), c_arr(nn, True), c_int_p, c_int_p, c_real_p] 1533 self.c_calc_comp.append(define_func(t, i, npl)) 1534 c_calc_pot_update = lib.calc_inv_update 1535 c_calc_pot_update.restype = None 1536 c_calc_pot_update.argtypes = [c_arr(npl)] 1537 c_get_dc = lib.get_dc 1538 c_get_dc.restype = None 1539 c_get_dc.argtypes = [c_arr(nn, True), c_arr(nx, True), c_arr(nni, True), c_arr(no, True), c_arr(ni, True)] 1540 def calc(u, x, v): 1541 u = np.array(u, dtype=self.dtp) 1542 x = np.array(x, dtype=self.dtp) 1543 v = np.array(v, dtype=self.dtp) 1544 x_new = np.zeros(nx, dtype=self.dtp) 1545 o = np.zeros(no, dtype=self.dtp) 1546 if c_calc(u, x, v, x_new, o, ct.byref(info), ct.byref(nfev), ct.byref(fnorm)) < 0: 1547 if fnorm.value > 1e-7: 1548 raise RuntimeError("convergence: method=%s, info=%d, nfev=%d, fnorm=%g" 1549 % (self.method, info.value, nfev.value, fnorm.value)) 1550 return o, x_new, v 1551 def calc_stream(u, ii=-1): 1552 assert u.shape[1] == ni + (ii >= 0) 1553 u = np.array(u, dtype=self.dtp, order="C", copy=False) 1554 o = np.zeros((u.shape[0], no), dtype=self.dtp) 1555 if c_calc_stream(u, o, u.shape[0], ii) != 0: 1556 v, x, minmax, g_info, g_nfev, g_fnorm = get_info() 1557 raise ValueError("convergence error: info=%d, nfev=%d, fnorm=%g" % (g_info, g_nfev, g_fnorm)) 1558 return o 1559 def get_info(): 1560 v = np.zeros(nn, dtype=self.dtp, order='C') 1561 x = np.zeros(nx, dtype=self.dtp, order='C') 1562 minv = np.zeros(nni, dtype=self.dtp, order='C') 1563 maxv = np.zeros(nni, dtype=self.dtp, order='C') 1564 c_get_info(v, x, minv, maxv, ct.byref(info), ct.byref(nfev), ct.byref(fnorm)) 1565 return v, x, np.vstack((minv, maxv)).T, info.value, nfev.value, fnorm.value 1566 def set_state(v, x): 1567 v = np.array(v, dtype=self.dtp) 1568 x = np.array(x, dtype=self.dtp) 1569 c_set_state(v, x) 1570 def get_dc(): 1571 v0 = np.zeros(nn, dtype=self.dtp, order='C') 1572 x0 = np.zeros(nx, dtype=self.dtp, order='C') 1573 p0 = np.zeros(nni, dtype=self.dtp, order='C') 1574 o0 = np.zeros(no, dtype=self.dtp, order='C') 1575 op = np.zeros(ni, dtype=self.dtp, order='C') 1576 c_get_dc(v0, x0, p0, o0, op) 1577 return v0, x0, ml.matrix(p0).T, o0, op 1578 self.c_calc = calc 1579 self.c_set_state = set_state 1580 self.c_get_info = get_info 1581 self.c_calc_stream = calc_stream 1582 self.c_calc_pot_update = c_calc_pot_update 1583 self.c_get_dc = get_dc 1584 1585 def __call__(self, v_in, ii=-1): 1586 return self.eval_c(v_in, ii) 1587 1588 1589class Parser(object): 1590 # Nl = incidence matrix denoting the nodes which potentials are controlling the currents 1591 # Nr = incidence matrix denoting the nodes where the current flows in (> 0) or out (< 0) 1592 # 1593 def __init__(self, S, V, fs, TR=True, create_filter=False, symbolic=False): 1594 self.fs = fs 1595 self.TR = TR # True: TR (trapezoidal) integration, False: BE backward euler 1596 self.create_filter = create_filter 1597 self.symbolic = symbolic 1598 self.mm = 2.0 if TR else 1.0 1599 self.update(S, V) 1600 1601 def update(self, S, V): 1602 self.V = V 1603 self.nodes = {} 1604 self.element_name = collections.defaultdict(list) 1605 tc = self.collect(S, V) 1606 n = len(self.nodes) + tc["V"] 1607 self.S = ml.zeros((n,n)) 1608 self.N = dict([(t, ml.zeros((tc[t[0]], n))) 1609 for t in "R","Xl","Xr","Nl","Nr","I","O","P"]) 1610 self.Pv = np.zeros(tc["P"]) 1611 self.pot_func = np.array((None,)*tc["P"]) 1612 self.ConstVoltages = ml.zeros(n) 1613 self.Z = np.zeros(tc["X"]) 1614 self.CZ = np.ones(tc["N"], dtype=int) 1615 self.f = np.array((None,)*tc["N"]) 1616 if self.create_filter or self.symbolic: 1617 alpha = sp.symbols("s") 1618 self.S = sp.Matrix(self.S) 1619 for k in self.N.keys(): 1620 self.N[k] = sp.Matrix(self.N[k]) 1621 if self.symbolic: 1622 self.Pv = sp.Matrix(self.Pv) 1623 d = {} 1624 for k, v in V.items(): 1625 if isinstance(k, Node): 1626 sym = sp.Symbol(str(k)) 1627 if isinstance(v, dict): 1628 v = v.copy() 1629 v["value"] = sym 1630 else: 1631 v = sym 1632 d[k] = v 1633 V = d 1634 else: 1635 alpha = self.mm * self.fs 1636 def map_conn(c): 1637 if isinstance(c, Out): 1638 return c 1639 if c is None: 1640 return None 1641 return self.nodes.get(c, -1) 1642 for row in S: 1643 row[0].process(self, [map_conn(c) for c in row[1:]], V.get(row[0]), alpha) 1644 self.pot = V.get("POT", {}) 1645 self.op = V.get("OP",[0.]*tc["I"]) 1646 self.pot_list = [] 1647 for a, f in self.pot_func: 1648 s = str(a) 1649 if s not in self.pot_list: 1650 self.pot_list.append(s) 1651 1652 def get_status(self): 1653 return "" 1654 1655 def matches(self, create_filter, symbolic): 1656 return self.create_filter == create_filter and self.symbolic == symbolic 1657 1658 def get_nonlin_funcs(self): 1659 return self.f 1660 1661 def set_function(self, idx, expr, vl, base): 1662 assert self.f[idx] is None 1663 self.f[idx] = (expr, vl, range(base, base+len(vl))) 1664 1665 def get_pot_funcs(self): 1666 return self.pot_func 1667 1668 def get_pot_attr(self): 1669 attrlist = [] 1670 for e in set([e[0] for e in self.element_name["P"]]): 1671 t = self.V[e] 1672 if not isinstance(t, dict): 1673 t = dict(value=t) 1674 var = t.get('var') 1675 if var is None: 1676 var = str(e)+"v" 1677 name = t.get('name', var) 1678 loga = t.get('a', 0) 1679 inv = t.get('inv', 0) 1680 if loga: 1681 if inv: 1682 expr = lambda a: (math.exp(loga * (1 - a)) - 1) / (math.exp(loga) - 1) 1683 else: 1684 expr = lambda a: (math.exp(loga * a) - 1) / (math.exp(loga) - 1) 1685 else: 1686 if inv: 1687 expr = lambda a: 1-a 1688 else: 1689 expr = lambda a: a 1690 attrlist.append((var, name, loga, inv, expr)) 1691 return attrlist 1692 1693 def get_variable_defaults(self): 1694 return dict([(a, self.pot.get(str(a), 0.5)) for a, f in self.pot_func]) 1695 1696 def extra_variable_index(self, idx): 1697 return len(self.nodes) + idx 1698 1699 def extra_variable_by_name(self, tpl): 1700 try: 1701 return len(self.nodes) + self.element_name["V"].index(tpl) 1702 except ValueError: 1703 print "%s not in %s" % (tpl, self.element_name["V"]) 1704 raise 1705 1706 def collect(self, S, V): 1707 tc = collections.Counter() 1708 l = self.element_name["S"] 1709 for row in S: 1710 e = row[0] 1711 conn = row[1:] 1712 e.add_count(tc, conn, V.get(e)) 1713 for s in conn: 1714 if s not in self.nodes: 1715 if s != GND and s is not None and not isinstance(s, Out): 1716 self.nodes[s] = len(self.nodes) 1717 l.append((s,None)) 1718 return tc 1719 1720 def new_row(self, N, sym, f=None): 1721 l = self.element_name[N] 1722 n = len(l) 1723 if N == "V": 1724 n += len(self.nodes) 1725 l.append((sym,f)) 1726 return n 1727 1728 def current_row(self, N): 1729 n = len(self.element_name[N]) - 1 1730 if N == "V": 1731 n += len(self.nodes) 1732 return n 1733 1734 def add_currents(self, m, conn, value): 1735 if conn[0] != -1: 1736 m[conn[0], conn[0]] += value 1737 if conn[1] != -1: 1738 m[conn[0], conn[1]] += -value 1739 if conn[1] != -1: 1740 m[conn[1], conn[1]] += value 1741 if conn[0] != -1: 1742 m[conn[1], conn[0]] += -value 1743 1744 def add_S_currents(self, conn, value): 1745 self.add_currents(self.S, conn, value) 1746 1747 def add_S(self, idx, conn, value): 1748 for i in conn[0], conn[1]: 1749 if i != -1: 1750 self.S[idx,i] += value 1751 value = -value 1752 for i in conn[0], conn[1]: 1753 value = -value 1754 if i != -1: 1755 self.S[i,idx] += value 1756 1757 def add_conn(self, N, node, conn, val): 1758 if conn != -1: 1759 self.N[N][node, conn] += val 1760 1761 def add_2conn(self, N, node, conn, value=1): 1762 if len(conn) != 2: 1763 raise ValueError("2 connections expected") 1764 for i, v in zip(conn, (value, -value)): 1765 self.add_conn(N, node, i, v) 1766 1767 @staticmethod 1768 def format_element(el, pref=""): 1769 n, d = el 1770 return "%s%s%s" % (pref, n, d and ("[%s]" % d) or "") 1771 1772 def node_names(self): 1773 return [self.format_element(v) for v in self.element_name["S"]+self.element_name["V"]] 1774 1775 def out_labels(self): 1776 return [str(v) for v in self.element_name["O"][0][0].conn] 1777 1778 def generate_c_code(self, d): 1779 # UNUSED 1780 def mk_sym(s, n): 1781 return sp.symbols(['%s[%d]' % (s, i) for i in range(len(self.element_name[n]))]) 1782 def ccode(d, ret, ex): 1783 r = ret + '[%d]' 1784 d[ret] = "\n".join([sp.ccode(ee, r % i) for i, ee in enumerate(ex)]).replace("\n","\n ") 1785 def ccodesum(d, ret, l): 1786 l2 = [m for m in l if m] 1787 if l2: 1788 ccode(d, ret, reduce(operator.add, l2)) 1789 else: 1790 d[ret] = '/* no %s */' % ret 1791 x = sp.Matrix(mk_sym('x', 'X')) 1792 u = sp.Matrix(mk_sym('u', 'I')) 1793 i = sp.Matrix(mk_sym('i', 'N')) 1794 v = mk_sym('v', 'N') 1795 p = mk_sym('p', 'N') 1796 # nonlinear function for root-finding 1797 l = [] 1798 for n, (expr, vl, base) in enumerate(self.eq.f): 1799 for j, e in enumerate(vl): 1800 expr = expr.subs(e, v[base+j]) 1801 l.append(expr) 1802 ccode(d, 'fvec', sp.Matrix(p) + sp.Matrix(self.K) * sp.Matrix(l) - sp.Matrix(np.diag(self.eq.CZ)) * sp.Matrix(v)) 1803 # "currents" 1804 ccode(d, 'i', l) 1805 # p value 1806 ccodesum(d, 'p', [sp.Matrix(self.G) * x, sp.Matrix(self.H) * u, sp.Matrix(self.Hc)]) 1807 # x update 1808 ccodesum(d, 'x_new', [sp.Matrix(self.A) * x, sp.Matrix(self.B) * u, sp.Matrix(self.Bc), sp.Matrix(self.C) * i]) 1809 # output 1810 ccodesum(d, 'o', [sp.Matrix(self.D) * x, sp.Matrix(self.E) * u, sp.Matrix(self.Ec), sp.Matrix(self.F) * i]) 1811 1812 1813def get_py_executor(parser, solver=None, linearize=False): 1814 sim = SimulatePy(EquationSystem(parser), solver) 1815 if linearize: 1816 J = sim.jacobi() 1817 sim = SimulatePy(EquationSystem(parser, J), solver) 1818 return sim 1819 1820def get_executor(name, parser, solver=None, pure_python=True, c_tempdir=None, c_verbose=False, 1821 c_debug_load="", c_real="double", extra_sources=None, linearize=False): 1822 if pure_python: 1823 return get_py_executor(parser, solver, linearize) 1824 elif c_debug_load: 1825 sim = SimulateC(c_debug_load) 1826 print "%s/%s: %s[%d], %s, %s" % (sim.name, sim.comment, sim.method, sim.data_size, sim.out_labels, sim.pot_list) 1827 return sim 1828 else: 1829 sim = SimulatePy(EquationSystem(parser), solver) 1830 return BuildCModule(name, sim, solver, c_tempdir, c_verbose, c_real, extra_sources, linearize).get_executor()[1] 1831