1 /*
2 * wavefile.c - stuff for working with entire datasets of waveform data.
3 *
4 * Copyright 1999, Stephen G. Tell.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public
17 * License along with this library; if not, write to the Free
18 * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
19 *
20 */
21
22 #include "ssintern.h"
23 #include <stdlib.h>
24 #include <stdio.h>
25 #include <string.h>
26 #include <errno.h>
27 #include <config.h>
28 #include <glib.h>
29 #include "wavefile.h"
30
31
32 #ifdef HAVE_POSIX_REGEXP
33 #include <regex.h>
34 #define REGEXP_T regex_t
35 #define regexp_test(c,s) (regexec((c), (s), 0, NULL, 0) == 0)
36 static regex_t *
regexp_compile(char * str)37 regexp_compile(char *str)
38 {
39 int err;
40 char ebuf[128];
41 regex_t *creg;
42
43 creg = g_new(regex_t, 1);
44 err = regcomp(creg, str, REG_NOSUB|REG_EXTENDED);
45 if(err) {
46 regerror(err, creg, ebuf, sizeof(ebuf));
47 fprintf(stderr, "internal error (in regexp %s):\n", str);
48 fprintf(stderr, " %s\n", ebuf);
49 exit(1);
50 }
51 return creg;
52 }
53
54 #else
55 #include "regexp.h" /* Henry Spencer's V8 regexp */
56 #define REGEXP_T regexp
57 #define regexp_test(c,s) regexec((c), (s))
58 #define regexp_compile(s) regcomp(s)
59 #endif
60
61 WaveFile *wf_finish_read(SpiceStream *ss);
62 WvTable *wf_read_table(SpiceStream *ss, WaveFile *wf, int *statep, double *ivalp, double *dvals);
63 void wf_init_dataset(WDataSet *ds);
64 inline void wf_set_point(WDataSet *ds, int n, double val);
65 void wf_free_dataset(WDataSet *ds);
66 WvTable *wvtable_new(WaveFile *wf);
67 void wt_free(WvTable *wt);
68
69 typedef struct {
70 char *name;
71 char *fnrexp;
72 REGEXP_T *creg;/* compiled form of regexp */
73 } DFormat;
74
75 /* table associating file typenames with filename regexps.
76 * Typenames should be those supported by spicefile.c.
77 *
78 * Filename patterns are full egrep-style
79 * regular expressions, NOT shell-style globs.
80 */
81 static DFormat format_tab[] = {
82 {"hspice", "\\.(tr|sw|ac)[0-9]$" },
83 {"cazm", "\\.[BNW]$" },
84 {"spice3raw", "\\.raw$" },
85 {"spice2raw", "\\.rawspice$" },
86 {"nsout", "\\.out$" },
87 {"ascii", "\\.(asc|acs|ascii)$" }, /* ascii / ACS format */
88 };
89 static const int NFormats = sizeof(format_tab)/sizeof(DFormat);
90
91 /*
92 * Read a waveform data file.
93 * If the format name is non-NULL, only tries reading in specified format.
94 * If format not specified, tries to guess based on filename, and if
95 * that fails, tries all of the readers until one sucedes.
96 * Returns NULL on failure after printing an error message.
97 *
98 * TODO: use some kind of callback or exception so that client
99 * can put the error messages in a GUI or somthing.
100 */
wf_read(char * name,char * format)101 WaveFile *wf_read(char *name, char *format)
102 {
103 FILE *fp;
104 SpiceStream *ss;
105 int i;
106
107 unsigned int tried = 0; /* bitmask of formats. */
108
109 g_assert(NFormats <= 8*sizeof(tried));
110 fp = fopen64(name, "r");
111 if(fp == NULL) {
112 perror(name);
113 return NULL;
114 }
115
116 if(format == NULL) {
117 for(i = 0; i < NFormats; i++) {
118 if(!format_tab[i].creg) {
119 format_tab[i].creg = regexp_compile(format_tab[i].fnrexp);
120 }
121 if(regexp_test(format_tab[i].creg, name))
122 {
123 tried |= 1<<i;
124 ss = ss_open_internal(fp, name, format_tab[i].name);
125 if(ss) {
126 ss_msg(INFO, "wf_read", "%s: read with format \"%s\"", name, format_tab[i].name);
127 return wf_finish_read(ss);
128 }
129
130 if(fseek(fp, 0L, SEEK_SET) < 0) {
131 perror(name);
132 return NULL;
133 }
134
135 }
136 }
137 if(tried == 0)
138 ss_msg(INFO, "wf_read", "%s: couldn't guess a format from filename suffix.", name);
139 /* no success with formats whose regexp matched filename,
140 * try the others.
141 */
142 for(i = 0; i < NFormats; i++) {
143 if((tried & (1<<i)) == 0) {
144 ss = ss_open_internal(fp, name, format_tab[i].name);
145 if(ss)
146 return wf_finish_read(ss);
147 tried |= 1<<i;
148 if(fseek(fp, 0L, SEEK_SET) < 0) {
149 perror(name);
150 return NULL;
151 }
152 }
153 }
154 ss_msg(ERR, "wf_read", "%s: couldn't read with any format\n", name);
155 return NULL;
156 } else { /* use specified format only */
157 ss = ss_open_internal(fp, name, format);
158 if(ss)
159 return wf_finish_read(ss);
160 else
161 return NULL;
162 }
163 }
164
165 /*
166 * read all of the data from a SpiceStream and store it in the WaveFile
167 * structure.
168 */
wf_finish_read(SpiceStream * ss)169 WaveFile *wf_finish_read(SpiceStream *ss)
170 {
171 WaveFile *wf;
172 double ival;
173 double *dvals;
174 WvTable *wt;
175 int state;
176 double *spar = NULL;
177
178 wf = g_new0(WaveFile, 1);
179 wf->ss = ss;
180 wf->tables = g_ptr_array_new();
181 dvals = g_new(double, ss->ncols);
182
183 state = 0;
184 do {
185 wt = wf_read_table(ss, wf, &state, &ival, dvals);
186 if(wt) {
187 ss_msg(DBG, "wf_finish_read", "table with %d rows; state=%d", wt->nvalues, state);
188 wt->swindex = wf->wf_ntables;
189 g_ptr_array_add(wf->tables, wt);
190 if(!wt->name) {
191 char tmp[128];
192 sprintf(tmp, "tbl%d", wf->wf_ntables);
193 wt->name = g_strdup(tmp);
194 }
195 } else {
196 ss_msg(DBG, "wf_finish_read", "NULL table; state=%d", state);
197 }
198 } while(state > 0);
199
200 g_free(dvals);
201 g_free(spar);
202 ss_close(ss);
203
204 if(state < 0) {
205 wf_free(wf);
206 return NULL;
207 } else {
208 return wf;
209 }
210 }
211
212 /*
213 * read data for a single table (sweep or segment) from spicestream.
214 * on entry:
215 * state=0: no previous data; dvals is allocated but garbage
216 * state=1: first row of data is in *ivalp, and vals[].
217 * on exit:
218 * return NULL: fatal error, *statep=-1
219 * return non-NULL: valid wvtable*
220 *
221 * state=-1 fatal error
222 * state=0: successful completion of reading whole file
223 * state=1: finished table but more tables remain,
224 * none of the next table has yet been read
225 * state=2: finished table but more tables remain and
226 * *ivalp,dvals[] contain first row of next table.
227 */
228 WvTable *
wf_read_table(SpiceStream * ss,WaveFile * wf,int * statep,double * ivalp,double * dvals)229 wf_read_table(SpiceStream *ss, WaveFile *wf,
230 int *statep, double *ivalp, double *dvals)
231 {
232 WvTable *wt;
233 int row;
234 WaveVar *dv;
235 double last_ival;
236 double spar;
237 int rc, i, j;
238
239 if(ss->nsweepparam > 0) {
240 if(ss->nsweepparam == 1) {
241 if(ss_readsweep(ss, &spar) <= 0) {
242 *statep = -1;
243 return NULL;
244 }
245 } else {
246 ss_msg(ERR, "wf_read_table", "nsweepparam=%d; multidimentional sweeps not supported\n", ss->nsweepparam);
247 *statep = -1;
248 return NULL;
249 }
250 }
251 wt = wvtable_new(wf);
252 if(ss->nsweepparam == 1) {
253 wt->swval = spar;
254 wt->name = g_strdup(ss->spar[0].name);
255 } else {
256 wt->swval = 0;
257 }
258
259 if(*statep == 2) {
260 wf_set_point(wt->iv->wds, row, *ivalp);
261 for(i = 0; i < wt->wt_ndv; i++) {
262 dv = &wt->dv[i];
263 for(j = 0; j < dv->wv_ncols; j++)
264 wf_set_point(&dv->wds[j], row,
265 dvals[dv->sv->col - 1 + j ]);
266 }
267 row = 1;
268 wt->nvalues = 1;
269 last_ival = *ivalp;
270 } else {
271 row = 0;
272 wt->nvalues = 0;
273 last_ival = -1.0e29;
274 }
275
276 while((rc = ss_readrow(ss, ivalp, dvals)) > 0) {
277 if(row > 0 && *ivalp < last_ival) {
278 if(row == 1) {
279 ss_msg(ERR, "wavefile_read", "independent variable is not nondecreasing at row %d; ival=%g last_ival=%g\n", row, *ivalp, last_ival);
280 wt_free(wt);
281 *statep = -1;
282 return NULL;
283
284 } else {
285 *statep = 2;
286 return wt;
287 }
288 }
289 last_ival = *ivalp;
290 wf_set_point(wt->iv->wds, row, *ivalp);
291 for(i = 0; i < wt->wt_ndv; i++) {
292 dv = &wt->dv[i];
293 for(j = 0; j < dv->wv_ncols; j++)
294 wf_set_point(&dv->wds[j], row,
295 dvals[dv->sv->col - 1 + j ]);
296 }
297 row++;
298 wt->nvalues++;
299 }
300 if(rc == -2)
301 *statep = 1;
302 else if(rc < 0) {
303 wt_free(wt);
304 *statep = -1;
305 return NULL;
306 } else {
307 *statep = 0;
308 }
309 return wt;
310 }
311
312
313 /*
314 * Free all memory used by a WaveFile
315 */
316 void
wf_free(WaveFile * wf)317 wf_free(WaveFile *wf)
318 {
319 int i;
320 WvTable *wt;
321 for(i = 0; i < wf->tables->len; i++) {
322 wt = wf_wtable(wf, i);
323 wt_free(wt);
324 }
325 g_ptr_array_free(wf->tables, 0);
326 ss_delete(wf->ss);
327 g_free(wf);
328 }
329
wt_free(WvTable * wt)330 void wt_free(WvTable *wt)
331 {
332 int i;
333 for(i = 0; i < wt->wt_ndv; i++)
334 wf_free_dataset(wt->dv[i].wds);
335 g_free(wt->dv);
336 wf_free_dataset(wt->iv->wds);
337 g_free(wt->iv);
338 if(wt->name)
339 g_free(wt->name);
340 g_free(wt);
341 }
342
343 /*
344 * create a new, empty WvTable for a WaveFile
345 */
346 WvTable *
wvtable_new(WaveFile * wf)347 wvtable_new(WaveFile *wf)
348 {
349 WvTable *wt;
350 SpiceStream *ss = wf->ss;
351 int i, j;
352
353 wt = g_new0(WvTable, 1);
354 wt->wf = wf;
355 wt->iv = g_new0(WaveVar, 1);
356 wt->iv->sv = ss->ivar;
357 wt->iv->wtable = wt;
358 wt->iv->wds = g_new0(WDataSet, 1);
359 wf_init_dataset(wt->iv->wds);
360
361 wt->dv = g_new0(WaveVar, wf->ss->ndv);
362 for(i = 0; i < wf->wf_ndv; i++) {
363 wt->dv[i].wtable = wt;
364 wt->dv[i].sv = &ss->dvar[i];
365 wt->dv[i].wds = g_new0(WDataSet, wt->dv[i].sv->ncols);
366 for(j = 0; j < wt->dv[i].sv->ncols; j++)
367 wf_init_dataset(&wt->dv[i].wds[j]);
368 }
369 return wt;
370 }
371
372
373 /*
374 * initialize common elements of WDataSet structure
375 */
376 void
wf_init_dataset(WDataSet * ds)377 wf_init_dataset(WDataSet *ds)
378 {
379 ds->min = G_MAXDOUBLE;
380 ds->max = -G_MAXDOUBLE;
381
382 ds->bpsize = DS_INBLKS;
383 ds->bptr = g_new0(double *, ds->bpsize);
384 ds->bptr[0] = g_new(double, DS_DBLKSIZE);
385 ds->bpused = 1;
386 ds->nreallocs = 0;
387 }
388
389 /*
390 * free up memory pointed to by a DataSet, but not the dataset itself.
391 */
392 void
wf_free_dataset(WDataSet * ds)393 wf_free_dataset(WDataSet *ds)
394 {
395 int i;
396 for(i = 0; i < ds->bpused; i++)
397 if(ds->bptr[i])
398 g_free(ds->bptr[i]);
399 g_free(ds->bptr);
400 g_free(ds);
401 }
402
403 /*
404 * Iterate over all WaveVars in all sweeps/segments in the WaveFile,
405 * calling the function for each one.
406 */
407 void
wf_foreach_wavevar(WaveFile * wf,GFunc func,gpointer * p)408 wf_foreach_wavevar(WaveFile *wf, GFunc func, gpointer *p)
409 {
410 WvTable *wt;
411 WaveVar *wv;
412 int i, j;
413
414 for(i = 0; i < wf->wf_ntables; i++) {
415 wt = wf_wtable(wf, i);
416 for(j = 0; j < wf->wf_ndv; j++) {
417 wv = &wt->dv[j];
418 (func)(wv, p);
419 }
420 }
421 }
422
423 /*
424 * expand dataset's storage to add one more block.
425 */
426 void
wf_expand_dset(WDataSet * ds)427 wf_expand_dset(WDataSet *ds)
428 {
429 if(ds->bpused >= ds->bpsize) {
430 ds->bpsize *= 2;
431 ds->bptr = g_realloc(ds->bptr, ds->bpsize * sizeof(double*));
432 ds->nreallocs++;
433 }
434 ds->bptr[ds->bpused++] = g_new(double, DS_DBLKSIZE);
435 }
436
437 /*
438 * set single value in dataset. Probably can be inlined.
439 */
440 void
wf_set_point(WDataSet * ds,int n,double val)441 wf_set_point(WDataSet *ds, int n, double val)
442 {
443 int blk, off;
444 blk = ds_blockno(n);
445 off = ds_offset(n);
446 while(blk >= ds->bpused)
447 wf_expand_dset(ds);
448
449 ds->bptr[blk][off] = val;
450 if(val < ds->min)
451 ds->min = val;
452 if(val > ds->max)
453 ds->max = val;
454 }
455
456 /*
457 * get single point from dataset. Probably can be inlined.
458 */
459 double
wds_get_point(WDataSet * ds,int n)460 wds_get_point(WDataSet *ds, int n)
461 {
462 int blk, off;
463 blk = ds_blockno(n);
464 off = ds_offset(n);
465 g_assert(blk <= ds->bpused);
466 g_assert(off < DS_DBLKSIZE);
467
468 return ds->bptr[blk][off];
469 }
470
471 /*
472 * Use a binary search to return the index of the point
473 * whose value is the largest not greater than ival.
474 * if ival is equal or greater than the max value of the
475 * independent variable, return the index of the last point.
476 *
477 * Only works on independent-variables, which we require to
478 * be nondecreasing and have only a single column.
479 *
480 * Further, if there are duplicate values, returns the highest index
481 * that has the same value.
482 */
483 int
wf_find_point(WaveVar * iv,double ival)484 wf_find_point(WaveVar *iv, double ival)
485 {
486 WDataSet *ds = iv->wds;
487 double cval;
488 int a, b;
489 int n = 0;
490
491 a = 0;
492 b = iv->wv_nvalues - 1;
493 if(ival >= ds->max)
494 return b;
495 while(a+1 < b) {
496 cval = wds_get_point(ds, (a+b)/2);
497 /* printf(" a=%d b=%d ival=%g cval=%g\n", a,b,ival,cval); */
498 if(ival < cval)
499 b = (a+b)/2;
500 else
501 a = (a+b)/2;
502
503
504 g_assert(n++ < 32); /* > 2 ** 32 points? must be a bug! */
505 }
506 return a;
507 }
508
509 /*
510 * return the value of the dependent variable dv at the point where
511 * its associated independent variable has the value ival.
512 *
513 * FIXME:tell
514 * make this fill in an array of dependent values,
515 * one for each column in the specified dependent variable.
516 * This will be better than making the client call us once for each column,
517 * because we'll only have to search for the independent value once.
518 * (quick hack until we need support for complex and other multicolumn vars:
519 * just return first column's value.)
520 */
521 double
wv_interp_value(WaveVar * dv,double ival)522 wv_interp_value(WaveVar *dv, double ival)
523 {
524 int li, ri; /* index of points to left and right of desired value */
525 double lx, rx; /* independent variable's value at li and ri */
526 double ly, ry; /* dependent variable's value at li and ri */
527 WaveVar *iv;
528
529 iv = dv->wv_iv;
530
531 li = wf_find_point(iv, ival);
532 ri = li + 1;
533 if(ri >= dv->wv_nvalues)
534 return wds_get_point(dv->wds, dv->wv_nvalues-1);
535
536 lx = wds_get_point(&iv->wds[0], li);
537 rx = wds_get_point(&iv->wds[0], ri);
538 /* g_assert(lx <= ival); */
539 if(li > 0 && lx > ival) {
540 fprintf(stderr, "wv_interp_value: assertion failed: lx <= ival for %s: ival=%g li=%d lx=%g\n", dv->wv_name, ival, li, lx);
541 }
542
543 ly = wds_get_point(&dv->wds[0], li);
544 ry = wds_get_point(&dv->wds[0], ri);
545
546 if(ival > rx) { /* no extrapolation allowed! */
547 return ry;
548 }
549
550 return ly + (ry - ly) * ((ival - lx)/(rx - lx));
551 }
552
553 /*
554 * Find a named variable, return pointer to WaveVar
555 */
556 WaveVar *
wf_find_variable(WaveFile * wf,char * varname,int swpno)557 wf_find_variable(WaveFile *wf, char *varname, int swpno)
558 {
559 int i;
560 WvTable *wt;
561 if(swpno >= wf->wf_ntables)
562 return NULL;
563
564 for(i = 0; i < wf->wf_ndv; i++) {
565 wt = wf_wtable(wf, swpno);
566 if(0==strcmp(wt->dv[i].wv_name, varname))
567 return &wt->dv[i];
568 }
569 return NULL;
570 }
571