1import os 2from libc.string cimport const_char, strcmp 3from libc.stdlib cimport malloc, free 4import numpy as np 5cimport numpy as np 6np.import_array() 7cdef extern from "time.h": 8 cdef struct timespec: 9 int tv_sec 10 int tv_nsec 11 ctypedef int clockid_t 12 int CLOCK_MONOTONIC 13 int clock_gettime(clockid_t clk_id, timespec *tp) 14 15cdef extern from "circuit.hpp": 16 ctypedef struct N_Vector: 17 pass 18 N_Vector N_VNew_Serial(long int vec_length) 19 void N_VDestroy_Serial(N_Vector v) 20 double *NV_DATA_S(N_Vector v) 21 int NV_LENGTH_S(N_Vector v) 22 23 int FS "fs" 24 25 cppclass UserData "ComponentBase::UserData": 26 double *inval 27 double *state 28 29 cppclass ComponentBase: 30 int NEQ 31 int NDIM 32 int NVALS 33 int N_IN 34 int N_OUT 35 int n_params 36 double *params 37 char **param_names 38 char **in_names 39 char **out_names 40 char **state_names 41 char **var_names 42 int *ix 43 int verbose 44 int func(N_Vector x, N_Vector u, UserData *user_data) 45 void update(N_Vector y, N_Vector x, N_Vector state) 46 int startvalue(N_Vector x, N_Vector s, N_Vector u) 47 int findzero(double *x, double *s, N_Vector u) 48 int set_range(int i, double lower, double upper, int size) 49 void init(N_Vector u0, char* fname) 50 int calc(N_Vector x, N_Vector s, N_Vector u) 51 52 cppclass TriodeCircuit(ComponentBase): 53 pass 54 55 cppclass CoupledTriodeCircuit(ComponentBase): 56 pass 57 58 cppclass PowerAmpGate(ComponentBase): 59 pass 60 61 cppclass PowerAmpPlate(ComponentBase): 62 pass 63 64 cppclass PhaseSplitter(ComponentBase): 65 void feedback_coeff(double *b0, double *b1, double *a1, double pres) 66 67 cdef cppclass CSpline "Spline": 68 CSpline(double *, int, double, double) 69 double eval(double) 70 71fs = FS 72 73 74class ConvergenceError(ValueError): 75 76 def __init__(self, err, x, s, tag=None): 77 self.error = err 78 self.x = x 79 self.s = s 80 self.tag = tag 81 t = "%s: " % tag if tag else "" 82 v1 = ", ".join((str(v) for v in x)) 83 if s is not None: 84 v2 = "/[%s]" % ", ".join((str(v) for v in s)) 85 else: 86 v2 = "" 87 ValueError.__init__( 88 self, "%sKINSol error %d at [%s]%s" % (t, err, v1, v2)) 89 90 def __reduce__(self): 91 return (ConvergenceError, (self.error, self.x, self.s, self.tag)) 92 93 94cdef int find_idx(const_char *k, const_char **p, int n): 95 cdef int i 96 for i in range(n): 97 if strcmp(p[i], k) == 0: 98 return i 99 return -1 100 101cdef to_list(const_char **p, int n): 102 cdef int i 103 cdef list l = [] 104 cdef char *t 105 for i in range(n): 106 t = <char*>(p[i]) 107 l.append(t) 108 return l 109 110cdef inline double ts_diff(timespec ts1, timespec ts2): 111 cdef double df = ts1.tv_sec - ts2.tv_sec 112 return df * 1e9 + (ts1.tv_nsec - ts2.tv_nsec) 113 114cdef class CalcBase: 115 cdef ComponentBase *pbase 116 cdef int *idx 117 cdef int *capt_idx 118 cdef int capt_len 119 cdef list capt_vars 120 cdef np.ndarray _captured 121 cdef N_Vector _state 122 cdef np.ndarray state_min 123 cdef np.ndarray state_max 124 cdef N_Vector inp 125 cdef N_Vector outp 126 cdef double time_per_sample 127 128 property NDIM: 129 def __get__(self): 130 return self.pbase.NDIM 131 property NEQ: 132 def __get__(self): 133 return self.pbase.NEQ 134 property NVALS: 135 def __get__(self): 136 return self.pbase.NVALS 137 property N_IN: 138 def __get__(self): 139 return self.pbase.N_IN 140 property N_OUT: 141 def __get__(self): 142 return self.pbase.N_OUT 143 property param_names: 144 def __get__(self): 145 return to_list(self.pbase.param_names, self.pbase.n_params) 146 property in_names: 147 def __get__(self): 148 return to_list(self.pbase.in_names, self.pbase.N_IN) 149 property out_names: 150 def __get__(self): 151 return to_list(self.pbase.out_names, self.pbase.N_OUT) 152 property state_names: 153 def __get__(self): 154 return to_list(self.pbase.state_names, self.pbase.NDIM-self.pbase.N_IN) 155 property var_names: 156 def __get__(self): 157 return to_list(self.pbase.var_names, self.pbase.NEQ) 158 property state: 159 def __get__(self): 160 cdef int i 161 return np.array([NV_DATA_S(self._state)[i] for i in range(self.pbase.NDIM-self.pbase.N_IN)]) 162 def __set__(self, v): 163 for i in range(self.pbase.NDIM-self.pbase.N_IN): 164 NV_DATA_S(self._state)[i] = v[i] 165 property min_state: 166 def __get__(self): 167 return self.state_min 168 property max_state: 169 def __get__(self): 170 return self.state_max 171 property nanosec_per_sample: 172 "time in nanoseconds measured for last compute call" 173 def __get__(self): 174 return self.time_per_sample 175 176 177 def __dealloc__(self): 178 free(self.idx) 179 free(self.capt_idx) 180 N_VDestroy_Serial(self.inp) 181 N_VDestroy_Serial(self._state) 182 N_VDestroy_Serial(self.outp) 183 del self.pbase 184 185 property capture_signals: 186 def __get__(self): 187 return self.capt_vars 188 def __set__(self, v): 189 cdef int i, j 190 cdef char *t 191 if isinstance(v, basestring): 192 v = [v] 193 l = [] 194 for n in v: 195 t = n 196 j = find_idx(t, self.pbase.var_names, self.pbase.NEQ) 197 if j < 0: 198 raise ValueError("not found: %s", n) 199 l.append(j) 200 i = len(l) 201 if i != self.capt_len: 202 free(self.capt_idx) 203 self.capt_idx = <int*>malloc(i * sizeof(int)) 204 self.capt_len = i 205 for j in range(i): 206 self.capt_idx[j] = l[j] 207 self.capt_vars = v 208 209 property captured: 210 def __get__(self): 211 return self._captured 212 213 cdef int set_var(self, char *k, double v): 214 for i in range(self.pbase.n_params): 215 if strcmp(k, self.pbase.param_names[i]) == 0: 216 self.pbase.params[i] = v 217 return 1 218 return 0 219 220 def init(self): 221 cdef int i 222 for k, v in self.circuit.items(): 223 if not self.set_var(k, v): 224 raise KeyError("unknown: %s" % k) 225 for i, l in enumerate(self.start_grid): 226 if not self.pbase.set_range(i, l[0], l[1], l[2]): 227 raise ValueError("can't set range nr %d" % i) 228 if len(self.u0) != self.pbase.NEQ: 229 raise ValueError("length %d expected for u0" % self.pbase.NEQ) 230 for i, v in enumerate(self.u0): 231 NV_DATA_S(self.outp)[i] = v 232 if "AMPSIM_CACHE" in os.environ: 233 fname = "%s.cache" % self.comp_id 234 else: 235 fname = None 236 if fname and os.path.exists(fname) and os.stat(fname).st_mtime < os.stat("circuit/components.py").st_mtime: 237 print "remove", fname 238 os.remove(fname) 239 cdef char *_fname = NULL 240 if fname: 241 _fname = fname 242 self.pbase.init(self.outp, _fname) 243 244 def pre_call(self, A, with_state): 245 return A 246 247 def post_call(self, O, with_state): 248 return O 249 250 cdef np.ndarray prepare_input(self, object a, int with_state, int* n_in, int* n_out, int* ndim): 251 cdef np.ndarray A = np.PyArray_FROMANY(a, np.NPY_DOUBLE, 0, 2, np.NPY_ALIGNED) 252 ndim[0] = np.PyArray_NDIM(A) 253 n_in[0] = self.pbase.NDIM if with_state else self.pbase.N_IN 254 n_out[0] = self.pbase.N_OUT 255 if with_state: 256 n_out[0] += self.pbase.NDIM-self.pbase.N_IN 257 if ndim[0] == 0: 258 if n_in[0] > 1: 259 raise ValueError("input array: bad shape") 260 A = np.PyArray_Reshape(A, (1, 1)) 261 if n_out[0] > 1: 262 ndim[0] = 1 263 elif ndim[0] == 1: 264 if n_in[0] == 1: 265 A = np.PyArray_Reshape(A, ((np.PyArray_DIM(A, 0), 1))) 266 if n_out[0] > 1: 267 ndim[0] = 2 268 elif np.PyArray_DIM(A, 0) == n_in[0]: 269 A = np.PyArray_Reshape(A, (1, n_in[0])) 270 if n_out[0] == 1: 271 ndim[0] = 0 272 else: 273 raise ValueError("input array: bad shape") 274 elif np.PyArray_NDIM(A) == 2: 275 if np.PyArray_DIM(A, 1) != n_in[0]: 276 raise ValueError("input array: bad shape %d/%d" % (np.PyArray_DIM(A, 1), n_in[0])) 277 if n_out[0] == 1: 278 ndim[0] = 1 279 return A 280 281 def __call__(self, a, int with_state=False): 282 cdef int ndim = 0 283 cdef int n_in = 0 284 cdef int n_out = 0 285 cdef np.ndarray A = self.prepare_input(a, with_state, &n_in, &n_out, &ndim) 286 cdef np.npy_intp dim[2] 287 dim[0] = np.PyArray_DIM(A, 0) 288 dim[1] = n_out 289 cdef np.ndarray O = np.PyArray_SimpleNew(2, dim, np.NPY_DOUBLE) 290 cdef np.flatiter itc 291 if self.capt_len: 292 dim[0] = np.PyArray_DIM(A, 0) 293 dim[1] = self.capt_len 294 self._captured = np.PyArray_SimpleNew(2, dim, np.NPY_DOUBLE) 295 itc = np.PyArray_IterNew(self._captured) 296 A = self.pre_call(A, with_state) 297 cdef int i, j, n 298 cdef np.flatiter it = np.PyArray_IterNew(A) 299 cdef np.flatiter ito = np.PyArray_IterNew(O) 300 cdef np.npy_intp ix[2] 301 cdef timespec t0, t1 302 cdef double v 303 clock_gettime(CLOCK_MONOTONIC, &t0) 304 for i in range(np.PyArray_DIM(A, 0)): 305 ix[0] = i 306 for j in range(self.pbase.N_IN): 307 ix[1] = j 308 np.PyArray_ITER_GOTO(it, ix) 309 NV_DATA_S(self.inp)[j] = (<double*>np.PyArray_ITER_DATA(it))[0] 310 if with_state: 311 for j in range(n_in - self.pbase.N_IN): 312 ix[1] = self.pbase.N_IN + j 313 np.PyArray_ITER_GOTO(it, ix) 314 NV_DATA_S(self._state)[j] = (<double*>np.PyArray_ITER_DATA(it))[0] 315 n = self.pbase.calc(self.inp, self._state, self.outp) 316 if n: 317 raise ConvergenceError(n, A[i], None if with_state else self.state, getattr(self,"comp_id",None)) 318 if self.capt_len: 319 for j in range(self.capt_len): 320 (<double*>np.PyArray_ITER_DATA(itc))[j] = NV_DATA_S(self.outp)[self.capt_idx[j]] 321 np.PyArray_ITER_NEXT(itc) 322 self.pbase.update(self.outp, self.inp, self._state) 323 for j in range(self.pbase.N_OUT): 324 (<double*>np.PyArray_ITER_DATA(ito))[0] = NV_DATA_S(self.outp)[self.idx[j]] 325 np.PyArray_ITER_NEXT(ito) 326 if with_state: 327 for j in range(n_out-self.pbase.N_OUT): 328 (<double*>np.PyArray_ITER_DATA(ito))[0] = NV_DATA_S(self._state)[j] 329 np.PyArray_ITER_NEXT(ito) 330 else: 331 for j in range(self.pbase.NDIM - self.pbase.N_IN): 332 v = NV_DATA_S(self._state)[j] 333 if v > (<double*>np.PyArray_DATA(self.state_max))[j]: 334 (<double*>np.PyArray_DATA(self.state_max))[j] = v 335 if v < (<double*>np.PyArray_DATA(self.state_min))[j]: 336 (<double*>np.PyArray_DATA(self.state_min))[j] = v 337 clock_gettime(CLOCK_MONOTONIC, &t1) 338 self.time_per_sample = ts_diff(t1,t0)/np.PyArray_DIM(A, 0) 339 O = self.post_call(O, with_state) 340 cdef np.PyArray_Dims dm 341 if ndim != np.PyArray_NDIM(O): 342 if ndim == 1: 343 if np.PyArray_NDIM(O) == 0: 344 dim[0] = 1 345 else: 346 dim[0] = np.PyArray_DIM(O, 0) * np.PyArray_DIM(O, 1) 347 else: 348 if np.PyArray_NDIM(O) == 0: 349 dim[0] = 1 350 dim[1] = 1 351 else: 352 dim[0] = np.PyArray_DIM(O, 0) 353 dim[1] = 1 354 dm.ptr = dim 355 dm.len = ndim 356 O = np.PyArray_Newshape(O, &dm, np.NPY_ANYORDER) 357 return O 358 359 def func(self, a, state): 360 cdef int i 361 cdef N_Vector x = N_VNew_Serial(self.pbase.N_IN) 362 cdef N_Vector s = N_VNew_Serial(self.pbase.NDIM-self.pbase.N_IN) 363 cdef N_Vector u = N_VNew_Serial(self.pbase.NEQ) 364 for i in range(self.pbase.N_IN): 365 NV_DATA_S(x)[i] = a[i] 366 for i in range(self.pbase.NDIM-self.pbase.N_IN): 367 NV_DATA_S(s)[i] = state[i] 368 i = self.pbase.calc(x, s, u) 369 if i: 370 raise ConvergenceError(i, a, state, getattr(self,"comp_id",None)) 371 self.pbase.update(u, x, s) 372 for i in range(self.pbase.NDIM-self.pbase.N_IN): 373 state[i] = NV_DATA_S(s)[i] 374 ret = [NV_DATA_S(u)[self.idx[i]] for i in range(self.pbase.N_OUT)] 375 N_VDestroy_Serial(x) 376 N_VDestroy_Serial(s) 377 N_VDestroy_Serial(u) 378 return ret 379 380 def startvalue(self, x, s): 381 cdef int i 382 cdef N_Vector X = N_VNew_Serial(self.pbase.N_IN) 383 cdef N_Vector S = N_VNew_Serial(self.pbase.NDIM-self.pbase.N_IN) 384 cdef N_Vector U = N_VNew_Serial(self.pbase.NEQ) 385 for i in range(self.pbase.N_IN): 386 NV_DATA_S(X)[i] = x[i] 387 for i in range(self.pbase.NDIM-self.pbase.N_IN): 388 NV_DATA_S(S)[i] = s[i] 389 self.pbase.startvalue(X, S, U) 390 ret = [NV_DATA_S(U)[i] for i in range(self.pbase.NEQ)] 391 N_VDestroy_Serial(X) 392 N_VDestroy_Serial(S) 393 N_VDestroy_Serial(U) 394 return ret 395 396 def findzero(self, x, s, u, verbose=False): 397 cdef int i 398 cdef N_Vector X = N_VNew_Serial(self.pbase.N_IN) 399 cdef N_Vector S = N_VNew_Serial(self.pbase.NDIM-self.pbase.N_IN) 400 cdef N_Vector U = N_VNew_Serial(self.pbase.NEQ) 401 for i in range(self.pbase.N_IN): 402 NV_DATA_S(X)[i] = x[i] 403 for i in range(self.pbase.NDIM-self.pbase.N_IN): 404 NV_DATA_S(S)[i] = s[i] 405 for i in range(self.pbase.NEQ): 406 NV_DATA_S(U)[i] = u[i] 407 self.pbase.verbose = verbose 408 i = self.pbase.findzero(NV_DATA_S(X), NV_DATA_S(S), U) 409 self.pbase.verbose = False 410 if i != 0: 411 raise ValueError 412 ret = [NV_DATA_S(U)[i] for i in range(self.pbase.NEQ)] 413 N_VDestroy_Serial(X) 414 N_VDestroy_Serial(S) 415 N_VDestroy_Serial(U) 416 return ret 417 418 def equations(self, x, s, u): 419 cdef UserData D 420 cdef N_Vector X = N_VNew_Serial(self.pbase.N_IN) 421 cdef N_Vector S = N_VNew_Serial(self.pbase.NDIM-self.pbase.N_IN) 422 cdef N_Vector U = N_VNew_Serial(self.pbase.NEQ) 423 cdef N_Vector F = N_VNew_Serial(self.pbase.NEQ) 424 for i in range(self.pbase.N_IN): 425 NV_DATA_S(X)[i] = x[i] 426 for i in range(self.pbase.NDIM-self.pbase.N_IN): 427 NV_DATA_S(S)[i] = s[i] 428 for i in range(self.pbase.NEQ): 429 NV_DATA_S(U)[i] = u[i] 430 D.inval = NV_DATA_S(X) 431 D.state = NV_DATA_S(S) 432 self.pbase.func(U, F, &D) 433 ret = [NV_DATA_S(F)[i] for i in range(self.pbase.NEQ)] 434 N_VDestroy_Serial(X) 435 N_VDestroy_Serial(S) 436 N_VDestroy_Serial(U) 437 N_VDestroy_Serial(F) 438 return ret 439 440 def __init__(self): 441 cdef np.npy_intp i 442 cdef int j 443 self.inp = N_VNew_Serial(self.pbase.N_IN) 444 self._state = N_VNew_Serial(self.pbase.NDIM-self.pbase.N_IN) 445 i = self.pbase.NDIM-self.pbase.N_IN 446 self.state_min = np.PyArray_SimpleNew(1, &i, np.NPY_DOUBLE) 447 self.state_min[:] = 1e99 448 self.state_max = np.PyArray_SimpleNew(1, &i, np.NPY_DOUBLE) 449 self.state_max[:] = -1e99 450 self.outp = N_VNew_Serial(self.pbase.NEQ) 451 self.idx = <int*>malloc(2*sizeof(int)) 452 for i in range(self.pbase.N_OUT): 453 j = find_idx(self.pbase.out_names[i], self.pbase.var_names, self.pbase.NEQ) 454 if j < 0: 455 raise ValueError("not found: %s", self.pbase.out_names[i]) 456 self.idx[i] = j 457 self.state = self.operating_point 458 459 460cdef class tcParams(CalcBase): 461 def __cinit__(self): 462 self.pbase = new TriodeCircuit() 463 464 465cdef class ctcParams(CalcBase): 466 def __cinit__(self): 467 self.pbase = new CoupledTriodeCircuit() 468 469 470cdef class pagParams(CalcBase): 471 def __cinit__(self): 472 self.pbase = new PowerAmpGate() 473 474cdef class papParams(CalcBase): 475 def __cinit__(self): 476 self.pbase = new PowerAmpPlate() 477 478 def Xpre_call(self, A, with_state): 479 return np.append(A[:,0], A[:,1] + A[:,0], axis=1) 480 481 def post_call(self, O, with_state): 482 if with_state: 483 return np.append(O[:,:1] - O[:,1:2], O[:,2:], axis=1) 484 else: 485 return O[:,0] - O[:,1] 486 487 def func(self, a, state): 488 # expects [Ug1, Ug1+Ug2] 489 # returns Ua1 - Ua2 490 cdef int i 491 cdef N_Vector x = N_VNew_Serial(self.pbase.N_IN) 492 cdef N_Vector s = N_VNew_Serial(self.pbase.NDIM-self.pbase.N_IN) 493 cdef N_Vector u = N_VNew_Serial(self.pbase.NEQ) 494 for i in range(self.pbase.N_IN): 495 NV_DATA_S(x)[i] = a[i] 496 for i in range(self.pbase.NDIM-self.pbase.N_IN): 497 NV_DATA_S(s)[i] = state[i] 498 #NV_DATA_S(x)[1] -= NV_DATA_S(x)[0] #added 499 i = self.pbase.calc(x, s, u) 500 if i: 501 raise ConvergenceError(i, a, state) 502 self.pbase.update(u, x, s) 503 for i in range(self.pbase.NDIM-self.pbase.N_IN): 504 state[i] = NV_DATA_S(s)[i] 505 ret = [NV_DATA_S(u)[3] - NV_DATA_S(u)[4]] #changed 506 N_VDestroy_Serial(x) 507 N_VDestroy_Serial(s) 508 N_VDestroy_Serial(u) 509 return ret 510 511 512cdef class psParams(CalcBase): 513 def __cinit__(self): 514 self.pbase = new PhaseSplitter() 515 516 def feedback_coeff(self, double pres): 517 cdef double b0, b1, a1 518 (<PhaseSplitter*>self.pbase).feedback_coeff(&b0, &b1, &a1, pres) 519 return b0, b1, a1 520 521 522cdef class Spline: 523 cdef CSpline *spl 524 cdef np.ndarray y 525 526 def __cinit__(self, np.ndarray[np.float_t,ndim=1] y, double xstart, double h): 527 self.y = y 528 self.spl = new CSpline(<double*>np.PyArray_DATA(y), len(y), xstart, h) 529 530 def __dealloc__(self): 531 del self.spl 532 533 def __call__(self, np.ndarray[np.float_t,ndim=1] x): 534 cdef int i 535 cdef double v 536 out = np.empty_like(x) 537 for i in range(len(x)): 538 out[i] = self.spl.eval(x[i]) 539 return out 540 541# Local Variables: 542# compile-command: "rm -f pycircuit.so; python setup.py build_ext --inplace" 543# End: 544