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