1from __future__ import division 2import matplotlib 3# matplotlib.use('Qt4Agg') # set when you prefer to use the Qt backend over the Tk one 4import sympy, math, shutil, sys, os, numpy, pylab, warnings, tempfile, logging, argparse 5from cStringIO import StringIO 6import dk_simulator, models, circ, mk_netlist, dk_lib, simu, signals, generate_code 7from signals import Signal 8from dk_lib import CircuitException, error 9import scipy.optimize as opt 10 11try: 12 get_ipython 13except NameError: 14 pass 15else: 16 def _circuit_exception_handler(self, etype, value, tb, tb_offset=None): 17 print "error:", value 18 get_ipython().set_custom_exc((CircuitException,), _circuit_exception_handler) 19 20 21def eng_str(x): 22 y = abs(x) 23 exponent = int(math.log10(y)) 24 engr_exponent = exponent - exponent % 3 25 y /= 10**engr_exponent 26 if engr_exponent: 27 engr_exponent = 'e' + str(engr_exponent) 28 else: 29 engr_exponent = '' 30 return "%s%g%s" % ('-' if x < 0 else '', y, engr_exponent) 31 32def check_keywords(fname, kw, *keywords): 33 for k in kw: 34 if k not in keywords: 35 raise CircuitException("%s: unknown keyword '%s'" % (fname, k)) 36 37class NonlinComponent(object): 38 39 def __init__(self): 40 self.max_jacobi_error = None 41 self.max_jacobi = None 42 self.ptp = None 43 self.nonlin_range = None 44 self.max_error = None 45 self.basegrid = None 46 47 48class DC_Values(object): 49 def __init__(self, v, x, p, o): 50 self.v = v 51 self.x = x 52 self.p = p 53 self.o = o 54 def show(self): 55 with dk_lib.printoptions(linewidth=200): 56 print "v =", self.v 57 print "x =", self.x 58 print "p =", self.p.A1 59 print "o =", self.o 60 61SIM_PY, SIM_C, SIM_TABLE = range(3) 62 63class Circuit(object): 64 65 have_plot = False 66 67 def __init__(self, copy_from=None, FS=96000): 68 self.tempdir = None 69 self.tempdir_keep = False 70 self.backward_euler = False 71 self.solver = None 72 self.build_verbose = False 73 self.module_id = None 74 self.module_name = None 75 self.module_desc = None 76 self.c_datatype = "double" 77 self.table_datatype = "double" #"float" 78 self.pre_filter = None 79 self.post_filter = None 80 self.code_generator = None 81 self.plugindef = None 82 self.build_script = None 83 self.resample = True 84 self.FS = FS 85 self.use_sim = SIM_C 86 self.sig = Signal() 87 self._clear() 88 self._rmtree = shutil.rmtree 89 self.subcircuits = {} 90 if copy_from is not None: 91 self.backward_euler = copy_from.backward_euler 92 self.solver = copy_from.solver 93 self.FS = copy_from.FS 94 self.S = copy_from.S 95 self.V = copy_from.V 96 97 # 98 # internal routines 99 # 100 101 def _clear(self): 102 self.max_homotopy_iter = None 103 self.dc_method = "A" 104 self.sys_reduce_tol = 0 105 self.solver_params = None 106 self.grid_error_limit = None 107 self.transform_opts = dk_simulator.TransformOptions() 108 self._clear_all() 109 110 def _clear_calculated(self): 111 self.eq = None 112 self.sim_py = None 113 self.sim_c = None 114 self.sim_table = None 115 self.sim_filter = None 116 self.dc_values = None 117 self.components = None 118 self.ptp = None 119 self.minmax = None 120 self.E = None 121 self.basegrid = None 122 self.knot_positions = None 123 self.table_source = None 124 self.input_amplitude = None 125 self.sensitivity = None 126 127 def _clear_all(self): 128 self.S = None 129 self.V = None 130 self.loaded_filename = None 131 self.parser = None 132 self._clear_calculated() 133 134 def __del__(self): 135 self._remove_tempdir(True) 136 137 def _remove_tempdir(self, ignore_errors=False): 138 if self.tempdir is None: 139 return 140 if self.tempdir_keep: 141 return 142 self._rmtree(self.tempdir, ignore_errors=ignore_errors) 143 144 def _ensure_parser(self, symbolic=False, for_filter=False): 145 if self.parser is None or not self.parser.matches(symbolic, for_filter): 146 if self.S is None: 147 raise CircuitException("no netlist defined") 148 if self.V is None and not symbolic: 149 raise CircuitException("no netlist values defined") 150 self.parser = dk_simulator.Parser(self.S, self.V, self.FS, not self.backward_euler, for_filter, symbolic) 151 152 def _ensure_eq(self): 153 if self.eq is None: 154 self._ensure_parser() 155 self.eq = dk_simulator.EquationSystem(self.parser, opts=self.transform_opts) 156 157 def _ensure_sim_py(self): 158 if self.sim_py is None: 159 self._ensure_eq() 160 try: 161 self.sim_py = dk_simulator.SimulatePy(self.eq, self.solver, self.dc_method) 162 if self.sys_reduce_tol >= 0: 163 self.sim_py.balance_realization(self.sys_reduce_tol) 164 if (isinstance(self.eq.nonlin, dk_simulator.PartitionedNonlinEquations) 165 and self.transform_opts.decompose): 166 self.eq.F, self.eq.C = self.eq.nonlin.decompose_blocks(self.eq.F, self.eq.C) 167 except dk_simulator.ConvergenceError as e: 168 raise CircuitException(e) 169 170 def _ensure_dc_values(self): 171 if self.dc_values is None: 172 self._ensure_sim_py() 173 self.dc_values = DC_Values(self.sim_py.v00, self.sim_py.x0, self.sim_py.p0, self.sim_py.o0) 174 175 def _ensure_sim_c(self): 176 if self.sim_c is None: 177 self._ensure_dc_values() 178 if self.tempdir is None: 179 self.tempdir = tempfile.mkdtemp("_dk_compile") 180 modc = dk_simulator.BuildCModule( 181 self._get_module_id(), self.sim_py, c_tempdir=self.tempdir, 182 c_verbose=self.build_verbose, c_real=self.c_datatype, pre_filter=self.pre_filter, post_filter=self.post_filter) 183 if self.max_homotopy_iter is not None: 184 modc.max_homotopy_iter = self.max_homotopy_iter 185 if self.solver_params is not None: 186 modc.solver_params = self.solver_params 187 self.c_module, self.sim_c = modc.get_executor() 188 189 def _check_input_amplitude(self): 190 if self.input_amplitude is None: 191 raise CircuitException("please define input_amplitude") 192 193 def calc_range(self, signal): 194 self._ensure_sim_c() 195 s = self.make_signal_vector(signal) 196 self.ptp = self.sim_c(s.input_signal).ptp() 197 pot_arr = numpy.array([[0., 1.]] * len(self.sim_c.pot_list)) 198 if pot_arr.size: 199 self.minmax = numpy.append(pot_arr, self.sim_c.minmax, axis=0) 200 else: 201 self.minmax = self.sim_c.minmax 202 logger = logging.getLogger("range") 203 for i, r in enumerate(self.sim_c.minmax): 204 logger.debug("nonlin range %d: %g .. %g" % (i, r[0], r[1])) 205 #print self.ptp 206 #print repr(self.minmax) 207 208 def _ensure_range(self): 209 if self.minmax is None: 210 self._check_input_amplitude() 211 self.calc_range(self.sig(self.input_amplitude*self.sig.sweep())) 212 213 def _get_op(self): 214 if self.sim_py: 215 return self.sim_py.op 216 if self.sim_c: 217 return self.sim_c.op 218 else: 219 self._ensure_parser() 220 return self.parser.op 221 222 def _estimate_sensitivity(self, i_ptp, ii): 223 self._check_input_amplitude() 224 s = self.input_amplitude * self.sig.sweep() 225 pert = i_ptp * 1e-7 * self.sig.sine(freq=1000) 226 p = self.sig(sympy.Tuple(s, pert)) 227 self.sim_c.reset() 228 ch = len(self._get_op()) + 1 229 signal = self.sig.generate(p, self.FS, channels=ch) 230 out_t = self.sim_c(signal.input_signal, ii) 231 max_pert = max(abs(signal.signal[:,-1])) 232 self.sim_c.reset() 233 signal = self.make_signal_vector(self.sig(s)) 234 out_c = self.sim_c(signal.input_signal) 235 df = out_t-out_c 236 err = 2 * max(abs(df)) / out_c.ptp() 237 return err / max_pert 238 239 def _calc_ptp(self, idx): 240 self._ensure_sim_c() 241 self._ensure_range() 242 def calc(sgn): 243 def ff(p): 244 return -sgn * self.sim_c.nonlin_c(p[numpy.newaxis,:])[0,idx] 245 return -sgn * ff(opt.brute(ff, self.minmax, finish=None)) 246 return calc(1) - calc(-1) 247 248 def _ensure_sensitivity(self): 249 if self.sensitivity is None: 250 self._ensure_sim_c() 251 self.sensitivity = [ 252 self._estimate_sensitivity(self._calc_ptp(i), i) if self.basegrid[i] is not None else 1 253 for i in range(self.sim_c.nno)] 254 logger = logging.getLogger("approx") 255 logger.info("nonlin function sensitivity: %s" 256 % ", ".join(["%g" % v for v in self.sensitivity])) 257 258 def _ensure_max_error(self): 259 if self.E is None: 260 if not self.sim_c.nno: 261 self.E = [] 262 return 263 self._ensure_sensitivity() 264 lim = self.grid_error_limit 265 if lim is None: 266 lim = 1e-2 267 self.E = [lim / self.sensitivity[i] for i in range(self.sim_c.nno)] 268 269 #self._ensure_eq() 270 #jacobi_estimate_error = 10#1e-1 271 #maxerr = 1e-4 272 #J, vals = simu.estimate_max_jacobi(self.sim_c.nonlin, self.minmax, jacobi_estimate_error,# self.sim_c.nno) 273 #J = J[:,J.shape[1]-J.shape[0]:] ##FIXME 274 #J = numpy.matrix(J) 275 #dv = numpy.amax(numpy.append(abs(self.eq.F), abs(self.eq.F)*J*abs(self.eq.G)*abs(self.eq.C), axis=0), axis=0).A1 276 #self.E = maxerr * self.ptp / numpy.where(dv == 0, 1e-20, dv) 277 #print "E =", self.E 278 279 #grd_shape = vals.grd.shape 280 #numpoints = numpy.product(grd_shape[1:]) 281 #grd = vals.grd.reshape(grd_shape[0], numpoints) 282 #fnc = vals.values.reshape(self.eq.nonlin.nno, numpoints) 283 #self.cov = numpy.cov(grd, fnc)[:len(grd),len(grd):] 284 #print self.cov 285 286 def _ensure_basegrid(self): 287 if self.basegrid is None: 288 self._ensure_sim_c() 289 self.basegrid = [[[None,('s',4)]] * (self.sim_c.npl + self.sim_c.nni)] * self.sim_c.nno 290 291 def _ensure_knot_positions(self): 292 #if self.knot_positions is None: 293 # xx 294 pass 295 296 def _check_basegrid(self, basegrid, ns, nni, npl, nno): 297 if nno != len(basegrid): 298 error("basegrid needs %d rows (found %d)" 299 % (nno, len(basegrid)), "approx") 300 for i, row in enumerate(basegrid): 301 if row is None: 302 continue 303 if nni+npl != len(row): 304 if npl: 305 t = "%d parameters + %d inputs" % (nni, npl) 306 else: 307 t = "%d" % nni 308 error("basegrid[%d] has %d elements for function %s, input dimension: %s " 309 % (i, len(row), ns, t), "approx") 310 311 def _ensure_table_source(self): 312 if self.table_source is None: 313 if self.sim_c.nno == 0: 314 self.table_source = {} 315 self.maptype = None 316 return 317 self._ensure_eq() 318 o = StringIO() 319 inst = StringIO() 320 h = StringIO() 321 npl = self.sim_c.npl 322 tables = {} 323 i0v = self.sim_c.nonlin_c(numpy.append(numpy.matrix([0.5]*npl),self.sim_c.p0.T,axis=1))[0] 324 if len(self.sim_c.comp_sz) > 1: 325 for i, ((v_slice, p_slice, i_slice), ns) in enumerate(zip(self.sim_c.comp_sz, self.sim_c.comp_namespace)): 326 nni = p_slice.stop - p_slice.start 327 nno = i_slice.stop - i_slice.start 328 self._check_basegrid(self.basegrid[i_slice], ns, nni, npl, nno) 329 slc = range(npl) + range(p_slice.start+npl, p_slice.stop+npl) 330 class Comp: 331 comp_id = ns 332 comp_name = ns 333 ranges = self.minmax[slc] 334 basegrid = numpy.array([(g, None, None, e, True) for g, e in zip(self.basegrid, self.E)])[i_slice] 335 NVALS = nno 336 N_IN = nni+npl 337 NDIM = nni+npl 338 i0 = i0v[i_slice] 339 @staticmethod 340 def __call__(v, with_state): 341 return self.sim_c.c_calc_comp[i](v) 342 ##FIXME: maptype 343 self.maptype, spl = simu.TableGenerator.write_files(Comp(), o, inst, h) 344 tables[ns] = spl 345 else: 346 ns = "nonlin" 347 self._check_basegrid(self.basegrid, ns, self.eq.nonlin.nni, npl, self.eq.nonlin.nno) 348 class Comp: 349 comp_id = ns 350 comp_name = self._get_module_id() 351 ranges = self.minmax 352 basegrid = [(g, None, None, e, False) for g, e in zip(self.basegrid, self.E)] 353 NVALS = self.eq.nonlin.nno 354 N_IN = self.eq.nonlin.nni+npl 355 NDIM = self.eq.nonlin.nni+npl 356 i0 = i0v 357 @staticmethod 358 def __call__(v, with_state): 359 return self.sim_c.nonlin(v) 360 self.maptype, spl = simu.TableGenerator.write_files(Comp(), o, inst, h) 361 tables[ns] = spl 362 intpp_inst = "\n".join(set(inst.getvalue().split("\n"))) 363 self.table_source = dict(data_c=o.getvalue(), data_h=h.getvalue(), intpp_inst=intpp_inst, tables=tables) 364 365 def _build_sim_table(self, dev_interface=True): 366 self._ensure_sim_c() 367 self.sim_c.reset() 368 self._ensure_range() 369 self._ensure_max_error() 370 self._ensure_basegrid() 371 self._ensure_knot_positions() 372 self._ensure_table_source() 373 name = self._get_module_id() 374 self._ensure_sim_py() 375 if self.solver is None: 376 solver = {} 377 else: 378 solver = self.solver.copy() 379 solver["method"] = "table" 380 modc = dk_simulator.BuildCModule( 381 name, self.sim_py, solver=solver, 382 extra_sources=self.table_source, c_tempdir="gencode", 383 c_verbose=self.build_verbose, c_real=self.table_datatype, 384 pre_filter=self.pre_filter, post_filter=self.post_filter, 385 generator=self.code_generator) 386 modc.dev_interface = dev_interface 387 modc.resample = self.resample 388 modc.build_script = self.build_script 389 if self.solver_params is not None: 390 modc.solver_params = self.solver_params 391 if self.plugindef: 392 modc.plugindef = self.plugindef 393 return modc 394 395 def _ensure_sim_table(self): 396 if self.sim_table is None: 397 self.table_module, self.sim_table = self._build_sim_table().get_executor() 398 399 def _get_sim(self): 400 if self.use_sim == SIM_C: 401 self._ensure_sim_c() 402 return self.sim_c 403 elif self.use_sim == SIM_TABLE: 404 self._ensure_sim_table() 405 return self.sim_table 406 else: 407 assert self.use_sim == SIM_PY 408 self._ensure_sim_py() 409 return self.sim_py 410 411 def _ensure_filter(self, symbolic): 412 self._ensure_parser(symbolic=symbolic, for_filter=True) 413 if len(self.parser.get_nonlin_funcs()) > 0: 414 if symbolic: 415 raise CircuitException("circuit is nonlinear: symbolic formula generation not supported") 416 p = dk_simulator.Parser(self.S, self.V, self.FS, not self.backward_euler) 417 sim = dk_simulator.SimulatePy(dk_simulator.EquationSystem(p), self.solver, self.dc_method) 418 J = sim.jacobi() 419 else: 420 J = None 421 self.sim_filter = dk_simulator.LinearFilter(self.parser, J) 422 423 def _check_netlist(self): 424 if self.S is None: 425 raise CircuitException("no netlist loaded") 426 427 def _get_module_id(self, module_id=None): 428 if module_id is not None: 429 return module_id 430 if self.module_id is not None: 431 return self.module_id 432 if self.loaded_filename is not None: 433 return os.path.splitext(os.path.basename(self.loaded_filename))[0] 434 return "sim_module" 435 436 def _nodes_from_names(self, elements): 437 d = {} 438 elset = set() 439 for e in self.parser.element_name["N"]: 440 e = e[0] 441 elset.add(e) 442 d[str(e)] = e 443 l = [] 444 for e in elements: 445 if e in d: 446 l.append(d[e]) 447 elif e in elset: 448 l.append(e) 449 else: 450 raise CircuitException("%s unknown. nonlinear circuit elements: %s" 451 % (e, ", ".join(d.keys()))) 452 return l 453 454 def _nonlin_function_list(self, elements): 455 elements = self._nodes_from_names(elements) 456 el = {} 457 if elements: 458 for i, e in enumerate(self.parser.element_name["N"]): 459 if e[0] in elements: 460 e = e[0] 461 if e not in el: 462 el[e] = [] 463 el[e].append(i) 464 else: 465 for i, e in enumerate(self.parser.element_name["N"]): 466 e = e[0] 467 if e not in el: 468 el[e] = [] 469 el[e].append(i) 470 return el 471 472 # 473 # user interface 474 # 475 476 def show_status(self): 477 if self.loaded_filename: 478 print "circuit loaded from: %s" % self.loaded_filename 479 if self.S: 480 print "circuit element count: %d" % len(self.S) 481 if self.module_id: 482 print "module id: %s" % self.module_id 483 if self.tempdir: 484 print "temp dir: %s" % self.tempdir 485 if self.tempdir_keep: 486 print "keep temp dir: %s" % self.tempdir_keep 487 if self.parser: 488 print "Parser: %s" % self.parser.get_status() 489 if self.eq: 490 print "EquationSystem: %s" % self.eq.get_status() 491 if self.sim_c: 492 print "C Executor: %s, %d" % (self.sim_c.soname, self.sim_c.nx) 493 494 def make_signal_vector(self, signal): 495 return self.sig.generate(signal, self.FS, self._get_op()) 496 497 def print_netlist(self, values=False): 498 self._check_netlist() 499 for row in self.S: 500 print "%s: %s" % (row[0], ", ".join([str(v) for v in row[1:]])) 501 if values: 502 print 503 for e, v in sorted((str(e), v) for e, v in self.V.items()): 504 if isinstance(v, float): 505 v = eng_str(v) 506 print "%s = %s" % (e, v) 507 508 def read_gschem(self, filename, defs=None): 509 self._clear_all() 510 v = dict(vars(models)) 511 v["Transistors"] = circ.Transistors 512 v["Diodes"] = circ.Diodes 513 v["Tubes"] = circ.Tubes 514 v["math"] = math 515 if defs: 516 v.update(defs) 517 exec mk_netlist.read_netlist(filename) in v 518 self.S = v["S"] 519 self.V = v["V"] 520 self.loaded_filename = filename 521 522 def read_netlist(self, filename): 523 self._clear_all() 524 v = vars(models) 525 v["Transistors"] = circ.Transistors 526 v["Diodes"] = circ.Diodes 527 v["Tubes"] = circ.Tubes 528 v["math"] = math 529 with open(filename) as f: 530 exec f in v 531 self.S = v["S"] 532 self.V = v["V"] 533 self.loaded_filename = filename 534 535 def set_netlist(self, S, V=None): 536 self._clear_all() 537 self.S = S 538 self.V = V 539 540 def set_tempdir(self, path): 541 self._remove_tempdir() 542 self.tempdir = path 543 544 def keep_tempdir(self, keep=True): 545 self.tempdir_keep = keep 546 self._remove_tempdir() 547 548 def set_use(self, use): 549 self.use_sim = use 550 551 def set_samplerate(self, fs): 552 self._clear_calculated() 553 self.parser = None 554 self.FS = fs 555 556 def set_float(self, use_float=True): 557 self.c_datatype = "float" if use_float else "double" 558 559 def set_solver(self, solver): 560 self.solver = generate_code.solver_set_defaults(solver) 561 562 def load_module(self, filename, clear_all=True): 563 sim_c = dk_simulator.SimulateC(filename) 564 if not clear_all and self.FS != sim_c.fs: 565 raise CircuitException("Samplerate mismatch: %d / %d" % (self.FS, sim_c.FS)) 566 if clear_all: 567 self._clear_all() 568 self.FS = sim_c.fs 569 self.sim_c = sim_c 570 571 def set_module_id(self, module_id): 572 self.module_id = module_id 573 574 def set_in_nets(self, *connections): 575 self._check_netlist() 576 self.S = list(self.S) 577 l = (models.IN,) + tuple(connections) 578 for i, row in enumerate(self.S): 579 if row[0] == models.IN: 580 self.S[i] = l 581 break 582 else: 583 self.S.append(l) 584 self._clear_calculated() 585 if self.parser is not None: 586 self.parser.update(self.S, self.V) 587 588 def set_out_nets(self, *connections): 589 self._check_netlist() 590 self.S = list(self.S) 591 l = (models.OUT,) + tuple(connections) 592 for i, row in enumerate(self.S): 593 if row[0] == models.OUT: 594 self.S[i] = l 595 break 596 else: 597 self.S.append(l) 598 self._clear_calculated() 599 if self.parser is not None: 600 self.parser.update(self.S, self.V) 601 602 def add_element(self, element, connections, value): 603 self._check_netlist() 604 self.S = list(self.S) + [(element,)+tuple(connections)] 605 self.V[element] = value 606 self._clear_calculated() 607 if self.parser is not None: 608 self.parser.update(self.S, self.V) 609 return 610 611 def remove_element(self, element): 612 self._check_netlist() 613 self.S = list(self.S) 614 for i, row in enumerate(self.S): 615 if str(row[0]) == element: 616 del self.S[i] 617 self._clear_calculated() 618 if self.parser is not None: 619 self.parser.update(self.S, self.V) 620 return 621 raise CircuitException("%s not found int netlist" % element) 622 623 def remove_connected(self, net, exclude=None): 624 self._check_netlist() 625 if exclude is None: 626 exclude = {} 627 self.S = list(self.S) 628 comp = {net} 629 found = False 630 last_size = 0 631 def excluded_row(row): 632 for c in row[1:]: 633 if c in exclude and str(row[0]) in exclude[c]: 634 return True 635 return False 636 while len(comp) > last_size: 637 last_size = len(comp) 638 dl =[] 639 for i, row in enumerate(self.S): 640 for c in row[1:]: 641 if c in comp and not excluded_row(row): 642 comp |= set([j for j in row[1:] if j != models.GND]) 643 dl.append(i) 644 found = True 645 break 646 for i in reversed(dl): 647 del self.S[i] 648 if not found: 649 raise CircuitException("net %s not found" % net) 650 self._clear_calculated() 651 if self.parser is not None: 652 self.parser.update(self.S, self.V) 653 654 def split_circuit(self, parts): 655 for block, (cuts, inputs, outputs) in parts.items(): 656 c = Circuit(copy_from=self) 657 for net, els in cuts.items(): 658 c.remove_connected(net, exclude=cuts) 659 if inputs: 660 c.set_in_nets(*[p[0] for p in inputs]) 661 c.inputs = [p[1] for p in inputs] 662 if outputs: 663 c.set_out_nets(*[p[0] for p in outputs]) 664 c.outputs = [p[1] for p in outputs] 665 self.subcircuits[block] = c 666 667 def join_net(self, net, target_net): 668 self._check_netlist() 669 l = [] 670 found = False 671 for i, row in enumerate(self.S): 672 for j, c in enumerate(row[1:]): 673 if c == net: 674 l.append((i, j, c)) 675 elif c == target_net: 676 found = True 677 if not found: 678 raise CircuitException("target_net %s not found" % target_net) 679 if not l: 680 raise CircuitException("net %s not found" % net) 681 S = list(self.S) 682 for i, j, c in l: 683 S[i] = list(S[i]) 684 S[i][j+1] = target_net 685 self.S = S 686 self._clear_calculated() 687 if self.parser is not None: 688 self.parser.update(self.S, self.V) 689 690 def print_func_names(self): 691 self._ensure_parser() 692 print ", ".join(["%s:%s" % e for e in self.parser.element_name["N"]]) 693 694 def get_dc_values(self): 695 self._ensure_dc_values() 696 return self.dc_values 697 698 def linearize(self, *elements, **kw): 699 check_keywords("linearize", kw, "keep_dc") 700 keep_dc = kw.get('keep_dc', len(elements) != 0) 701 self._ensure_parser() 702 el = self._nonlin_function_list(elements) 703 S = [tuple([models.NODES]+[v[1] for v in sorted([(v, k) for k, v in self.parser.nodes.items()])])] 704 V = dict(self.V) 705 for e in self.S: 706 if e[0] in el: 707 if isinstance(e[0], models.CC_N): 708 e = (models.CC_L(e[0].d),) + e[1:] 709 elif isinstance(e[0], (models.Trans_F, models.Trans_GC)): 710 #FIXME calculate R (Trans_L parameter) 711 e = (models.Trans_L(e[0].d),) + e[1:] 712 elif isinstance(e[0], models.OPA): 713 p = V[e[0]] 714 if isinstance(p, dict): 715 p = p["A"] 716 V[e[0]] = p 717 S.append(e) 718 self.S = S 719 self.parser.update(S, V) 720 self.eq = None 721 self.sim_py = None 722 self.dc_values = None 723 self._ensure_dc_values() 724 el = self._nonlin_function_list(elements) 725 J = self.sim_py.jacobi() 726 Jc = self.sim_py.calc_i(self.dc_values.v) 727 S = [] 728 Nr = self.parser.N["Nr"] 729 Nl = self.parser.N["Nl"] 730 nodes = list(sorted(self.parser.nodes.keys(), key=lambda v: self.parser.nodes[v])) 731 def idx(a, v): 732 a = numpy.nonzero(numpy.ravel(a == v))[0] 733 if len(a) == 0: 734 return models.GND 735 else: 736 return nodes[int(a)] 737 for e in self.S: 738 if e[0] in el: 739 for i in el[e[0]]: 740 jl = numpy.ravel(J[i].nonzero()[1]) 741 for j in jl: 742 v = models.VCCS(str(e) + "l%d%d" % (i, j)) 743 S.append((v, idx(Nl[j], 1), idx(Nl[j], -1), idx(Nr[i], 1), idx(Nr[i], -1))) 744 dG = -J[i,j] 745 if keep_dc: 746 i0 = Jc[i]/len(jl) + dG*self.dc_values.v[j] 747 else: 748 i0 = 0 749 V[v] = dict(dG = dG, i0 = i0) 750 elif not keep_dc and isinstance(e[0], models.V): 751 V[e[0]] = 0 752 S.append(e) 753 else: 754 S.append(e) 755 self.S = S 756 self.V = V 757 self.parser.update(S, V) 758 self.eq = None 759 self.sim_py = None 760 self.dc_values = None 761 762 def print_filter_coeffs(self, symbolic=False, filename=None): 763 self._ensure_filter(symbolic=symbolic) 764 if filename is not None: 765 f = open(filename, "w") 766 else: 767 f = sys.stdout 768 b, a, terms = self.sim_filter.get_s_coeffs() 769 self.sim_filter.print_coeffs('b', b, f) 770 self.sim_filter.print_coeffs('a', a, f) 771 B, A, c = self.sim_filter.transform_bilinear(terms) 772 print >>f, "\nc = %s;" % c 773 self.sim_filter.print_coeffs('B', B, f) 774 self.sim_filter.print_coeffs('A', A, f) 775 if filename is not None: 776 f.close() 777 778 def get_faust_code(self, module_id=None, symbolic=False, filename=None, FS=None, pre_filter=None): 779 self._ensure_filter(symbolic=symbolic) 780 b, a = self.sim_filter.get_z_coeffs(samplerate=FS) 781 l = self.parser.get_pot_attr() 782 m_id = self._get_module_id(module_id) 783 plugindef = self.plugindef 784 if not plugindef: 785 plugindef = PluginDef(m_id) 786 dsp, ui = simu.generate_faust_module(plugindef, b, a, l, self.sim_filter, pre_filter, build_script=self.build_script) 787 return dsp, ui 788 789 def get_simple_faust_code(self, module_id=None, symbolic=False, filename=None, FS=None, pre_filter=None): 790 self._ensure_filter(symbolic=symbolic) 791 b, a = self.sim_filter.get_z_coeffs(samplerate=FS) 792 l = self.parser.get_pot_attr() 793 m_id = self._get_module_id(module_id) 794 plugindef = self.plugindef 795 if not plugindef: 796 plugindef = PluginDef(m_id) 797 dsp, ui = simu.generate_simple_faust_module(plugindef, b, a, l, self.sim_filter, pre_filter, build_script=self.build_script) 798 return dsp, ui 799 800 def save_faust_code(self, module_id=None, symbolic=False, filename=None, FS=None, pre_filter=None): 801 self._ensure_filter(symbolic=symbolic) 802 b, a = self.sim_filter.get_z_coeffs(samplerate=FS) 803 l = self.parser.get_pot_attr() 804 m_id = self._get_module_id(module_id) 805 plugindef = self.plugindef 806 if not plugindef: 807 plugindef = PluginDef(m_id) 808 dsp, ui = simu.generate_faust_module(plugindef, b, a, l, self.sim_filter, pre_filter, build_script=self.build_script) 809 if filename is None: 810 sys.stdout.write(dsp) 811 sys.stdout.write(ui) 812 else: 813 dspname = "{0}.dsp".format(filename) 814 uiname = "{0}_ui.cc".format(filename) 815 with open(dspname,"w") as f: 816 f.write(dsp) 817 with open(uiname,"w") as f: 818 f.write(ui) 819 820 def create_faust_module(self, module_id=None, symbolic=False, FS=None, pre_filter=None): 821 self._ensure_filter(symbolic=symbolic) 822 b, a = self.sim_filter.get_z_coeffs(samplerate=FS) 823 l = self.parser.get_pot_attr() 824 m_id = self._get_module_id(module_id) 825 plugindef = self.plugindef 826 if not plugindef: 827 plugindef = PluginDef(m_id) 828 simu.build_faust_module(plugindef, b, a, l, self.sim_filter, self.c_datatype, pre_filter, build_script=self.build_script) 829 830 def stream(self, signal, ii=-1): 831 if not isinstance(signal, signals.GeneratedSignal): 832 signal = self.make_signal_vector(signal) 833 sim = self._get_sim() 834 self.last_signal = signal 835 self.last_output = sim(signal.input_signal, ii) 836 return self.last_output 837 838 def plot(self, sig=None, label=None, clip=-80, nharmonics=8, spectrum=None, freq_range=None): 839 if sig is not None: 840 self.stream(sig) 841 if self.last_output is None: 842 raise CircuitException("nothing to plot") 843 lines = self.last_signal.plot(self.last_output, sig, label, clip, nharmonics, spectrum, freq_range) 844 #if label is None: 845 # label = self._get_sim().out_labels 846 #elif not hasattr(label, "__iter__"): 847 # l = [] 848 # for lbl in self._get_sim().out_labels: 849 # l.append("%s.%s" % (label, lbl)) 850 # label = l 851 Circuit.have_plot = True 852 return lines 853 854 def build_guitarix_module(self, include_paths=()): 855 modc = self._build_sim_table(dev_interface=False) 856 cw = modc.compile_c_func() 857 cw.script_dict["includes"] += "".join([" -I%s" % os.path.abspath(p) for p in include_paths]) 858 cw.wrap_code(load_module=False) 859 860 def deploy(self, path=None): 861 sim = self._get_sim() 862 if sim is self.sim_py: 863 self._ensure_sim_c() 864 sim = self.sim_c 865 if is_test(): 866 return 867 if sim is self.sim_table: 868 mod = self.table_module 869 s = "table" 870 else: 871 mod = self.c_module 872 s = "C" 873 if self.plugindef.lv2_plugin_type: 874 def try_mkdir(path): 875 try: 876 os.mkdir(path) 877 except OSError as e: 878 if e.errno != 17: # EEXIST 879 raise 880 vid = self.plugindef.lv2_versioned_id 881 if path is None: 882 path = os.path.expanduser("~/.lv2") 883 try_mkdir(path) 884 path = os.path.join(path, "gx_%s.lv2" % vid) 885 try_mkdir(path) 886 base = os.path.dirname(mod) 887 ttlname = os.path.join(base, vid+".ttl") 888 shutil.copy(ttlname, path) 889 shutil.copy(os.path.join(base, "manifest.ttl"), path) 890 fname = os.path.join(path, vid+".so") 891 s = "LV2 bundle [%s]" % s 892 shutil.copy(mod, fname) 893 print "LV2 bundle [%s module] copied to '%s'" % (s, path) 894 else: 895 if path is None: 896 path = os.path.expanduser("~/.config/guitarix/plugins/.") 897 fname = os.path.join(path, self._get_module_id()+".so") 898 shutil.copy(mod, fname) 899 print "%s module copied to '%s'" % (s, fname) 900 901 @staticmethod 902 def _get_samples(data, count): 903 return data[numpy.array(numpy.linspace(0, len(data)-1, count).round(), dtype=int)] 904 905 def check_result(self, sig, result, max_error=1e-7, count=None): 906 if not is_test(): 907 return True 908 test = get_test() 909 if count is None: 910 if result is None: 911 count = 10 912 else: 913 count = len(result) 914 y = self.stream(sig) 915 samples = self._get_samples(y, count) 916 if test.plot: 917 timeline = self.last_signal.timeline 918 pylab.plot(timeline, y) 919 if result is not None: 920 pylab.plot(self._get_samples(timeline, count), result, "rx") 921 pylab.show() 922 if test.printout: 923 print repr(samples) 924 return True 925 error = numpy.max(abs(result - samples)) / numpy.max(abs(result)) 926 if (error > max_error).any(): 927 print "%s: Difference = %g (> %g)" % (self._get_module_id(), error, max_error) 928 return False 929 else: 930 print "%s: OK" % self._get_module_id() 931 return True 932 933 def set_pot_variable(self, name, val): 934 self._get_sim().set_variable(name, val) 935 936 def set_pot_position(self, name, val): 937 self._ensure_parser() 938 expr = dict([(row[0], row[4]) for row in self.parser.get_pot_attr()])[name] 939 self.set_pot_variable(name, expr(val)) 940 941 def calc_nonlin_range(self, signal): 942 pass 943 944 def set_nonlin_range(self, idx, low=None, high=None): 945 pass 946 947 def estimate_max_jacobi(self, max_error): 948 pass 949 950 def set_max_jacobi(self, jacobi): 951 pass 952 953 def generate_table(self): 954 pass 955 956 957def _create_funcs(): 958 gcircuit = Circuit() 959 g = globals() 960 g['get_global_circuit'] = lambda: gcircuit 961 for v in dir(Circuit): 962 if not v.startswith("_"): 963 g[v] = getattr(gcircuit, v) 964_create_funcs() 965 966def show_plots(loc=None): 967 if Circuit.have_plot: 968 pylab.grid(which="both") 969 pylab.legend(loc=loc) 970 pylab.show() 971 Circuit.have_plot = False 972 973INFO = logging.INFO 974DEBUG = logging.DEBUG 975 976def set_log_level(lvl, logger=None): 977 logging.getLogger(logger).setLevel(lvl) 978 979class TestSettings: 980 active = False 981 plot = False 982 printout = False 983 984_display_traceback = True #False 985_testing = TestSettings 986 987def is_test(): 988 return _testing.active 989 990def get_test(): 991 return _testing 992 993def display_traceback(v): 994 global _display_traceback 995 _display_traceback = v 996 997def _init(): 998 global _testing 999 def _catch_circuit_error(tp, value, traceback): 1000 if tp == CircuitException: 1001 if _display_traceback: 1002 logging.getLogger(value.logger).error(value) 1003 else: 1004 print "ERROR:%s:%s" % (value.logger or "root", value) 1005 return 1006 if _excepthook is not None: 1007 _excepthook(tp, value, traceback) 1008 1009 logging.basicConfig() 1010 _excepthook = sys.excepthook 1011 sys.excepthook = _catch_circuit_error 1012 parser = argparse.ArgumentParser() 1013 parser.add_argument('--test', action='store_true', help='check results of tests') 1014 parser.add_argument('--test-plot', action='store_true', help='plot the test-data with markers') 1015 parser.add_argument('--test-print', action='store_true', help='print the result array') 1016 args = parser.parse_args() 1017 _testing.active = args.test or args.test_plot or args.test_print 1018 _testing.plot = args.test_plot 1019 _testing.printout = args.test_print 1020 1021_init() 1022