File: utils.i

package info (click to toggle)
yorick-yutils 1.3.0-2
  • links: PTS
  • area: main
  • in suites: lenny
  • size: 604 kB
  • ctags: 17
  • sloc: makefile: 105; python: 12
file content (1323 lines) | stat: -rw-r--r-- 47,762 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
/*
 * utils.i --
 *
 *	General purpose utility routines for Yorick.
 *
 * Copyright (c) 1996-2004, Eric Thibaut, Observatoire de Lyon (France).
 *
 * This program is free software; you can redistribute it and/or  modify it
 * under the terms of the GNU General Public License  as  published  by the
 * Free Software Foundation; either version 2 of the License,  or  (at your
 * option) any later version.
 *
 * This program is distributed in the hope  that  it  will  be  useful, but
 * WITHOUT  ANY   WARRANTY;   without   even   the   implied   warranty  of
 * MERCHANTABILITY or  FITNESS  FOR  A  PARTICULAR  PURPOSE.   See  the GNU
 * General Public License for more details (to receive a  copy  of  the GNU
 * General Public License, write to the Free Software Foundation, Inc., 675
 * Mass Ave, Cambridge, MA 02139, USA).
 * 
 * Routines:
 *      eval         - evaluate textual code.
 *	pwd          - print/get working directory.
 *	unref        - returns X, destroying X in the process.
 *	swap         - exchanges contents of two variables.
 *	grow_dimlist - build dimension list.
 *	reform       - change dimension list of an array.
 *	undersample  - shrink array dimension(s).
 *	spline_zoom  - resize an array by cubic spline interpolation
 *	map          - apply function to every elements of array or list.
 *	str(to)upper - convert string(s) to upper case letters.
 *	str(to)lower - convert string(s) to lower case letters.
 *	strcut       - cut text to fit into lines of given length.
 *	strjoin      - make array of strings into a single string.
 *	is_scalar    - check if argument is a scalar.
 *	is_vector    - check if argument is a vector.
 *	is_matrix    - check if argument is a matrix.
 *	is_real      - check if argument is of float or double type.
 *	is_complex   - check if argument is of complex type.
 *	is_integer   - check if argument is of integer type.
 *	is_numerical - check if argument is of numerical type.
 *	is_integer_scalar - check if argument is a scalar of integer type.
 *      expand_file_name - expand leading tilde in file name(s).
 *	get_file_name - get path name of file associated with a stream.
 *	protect_file_name - protect special characters in file name.
 *	read_ascii   - read ascii numeric data
 *	load_text    - load all lines of a text file.
 *	dump_text    - dump string array in a text file.
 *	guess_compression - guess compression method of existing file
 *	xopen        - extended open with (de)compression, primitives, ...
 *	stat         - display statistics/info about symbols/expressions.
 *	smooth       - smooth array along its dimensions.
 *      pdb_list     - lists contents of PDB binary file.
 *      pdb_restore_all - restore all non-record variables of a PDB file.
 *      timer_start  - (re)start the profiling timer.
 *      timer_elapsed - get/print the elapsed time since timer_start.
 *
 * History:
 *	$Id: utils.i,v 1.1.1.1 2007/12/11 23:55:13 frigaut Exp $
 *	$Log: utils.i,v $
 *	Revision 1.1.1.1  2007/12/11 23:55:13  frigaut
 *	Initial Import - yorick-yutils
 *	
 *	Revision 1.12  2004/10/14 09:51:37  eric
 *	 - Unused (to my knowledge) function "filenameof" removed,
 *	   it is superseded by "get_file_name".
 *	 - New functions: "expand_file_name", "get_file_name",
 *	   and "protect_file_name".
 *
 *	Revision 1.11  2004/08/31 07:26:17  eric
 *	 - New routines for PDB files: pdb_list, pdb_restore_all.
 *	 - New routines for profiling: timer_start, timer_elapsed.
 *
 *	Revision 1.10  2003/06/03 20:37:36  eric
 *	 - Function map() now accept a list argument (as _map) and use
 *	   weird names for local variables to avoid name clash.
 *	 - New functions: xopen, guess_compression, read_ascii,
 *	   load_text, and dump_text.
 *
 *	Revision 1.9  2002/11/22 08:24:27  eric
 *	 - slight changes in eval() code
 *	 - fix list of routines in leading comments of this file
 *
 *	Revision 1.8  2002/11/20 09:22:15  eric
 *	 - take care of not overwriting builtin functions (unref, strlower,
 *	   strupper, ...)
 *	 - new function: grow_dimlist
 *
 *	Revision 1.7  2002/07/25 10:52:32  eric
 *	- New function: eval.
 *
 *	Revision 1.6  2002/07/01 09:41:45  eric
 *	 - New function: pwd.
 *
 *	Revision 1.5  2002/06/06 14:19:42  eric
 *	 - New function: undersample.
 *
 *	Revision 1.4  2002/02/22 16:19:24  eric
 *	 - Change (one more time) names of str(to)upper/lower functions to
 *         avoid clash with builtin Yeti routines.
 *	 - Add "unref" routines (after Yorick's FAQ, so the true author is
 *         Dave, not me).
 *
 *	Revision 1.3  2001/12/08 22:33:37  eric
 *	 - stat(): computation of standard deviation improved.
 *
 *	Revision 1.2  2001/11/26 08:17:11  eric
 *	 - Functions strto{lower,upper} renamed as str{lower,upper}.
 *
 *	Revision 1.1  2001/03/23 16:20:52  eric
 *	Initial revision
 *
 *-----------------------------------------------------------------------------
 */

func eval(code, tmp=, debug=)
/* DOCUMENT eval, code;
       -or- eval(code);       
     Evaluates CODE given as a string or as an array of strings (considered
     as  different lines  in the  script).  Since  CODE can  be dynamically
     build,   this  routine   allows  the   execution  of   virtually  (see
     hints/restrictions below)  any Yorick's code  (e.g. dynamic definition
     of  structures,  of functions,  etc.).   For  instance, the  following
     statement defines a new structure:
       eval, "struct NewStruct {string a; long b; float c, d;}";

     Since  the script  gets evaluated  at the  scope level  of  the "eval"
     routine some local variables of the  "eval" routine may be used in the
     script:
       "eval_tmp"    contains  the  name of  the temporary script  file and
                     must not be changed by the script;
       "eval_debug"  contains the value of  the keyword DEBUG and  must not
                     be changed by the script;
       "eval_code"   contains the value of the argument CODE;
       "eval_result" is  returned by  "eval", its  contents may  be defined
                     into the script to provide a returned value.
     Note: impredictible  results may  occur if CODE  changes the  value of
     symbols "eval_tmp" and "eval_debug".

     Keyword TMP  can be  used to  specify the file  name of  the temporary
     script.  The default file name is:
       "$YORICK_EVAL_TMP"      if environment variable "YORICK_EVAL_TMP" is
                               set;
       "/tmp/$USER-eval_tmp.i" if environment variable "USER" set;
       "~/.eval_tmp.i"         otherwise.
     If  keyword DEBUG  is true  (non-zero and  non-nil), the  name  of the
     temporary file is printed out and the file is not removed.


   SEE ALSO: include. */
{
  /* Dump script into a temporary file. */
  if (is_void(tmp)) {
    /* Create default name for yorick temporary code. */
    tmp = get_env("YORICK_EVAL_TMP");
    if (! tmp) {
      tmp = get_env("USER");
      tmp = (tmp ? "/tmp/"+tmp+"-" : "~/.") + "eval_tmp.i";
    }
  }
  write, format="%s\n", open(tmp, "w"), code;
  
  /* Use "eval_" prefix in order to somewhat protect local variables
     from caller's code. */
  local eval_result, eval_tmp, eval_code, eval_debug;
  eq_nocopy, eval_tmp,   tmp;
  eq_nocopy, eval_debug, debug;
  eq_nocopy, eval_code,  code;

  /* Source script and return result. */
  include, eval_tmp, 1;
  if (eval_debug) write, format="Yorick code written in \"%s\"\n", eval_tmp;
  else remove, eval_tmp;
  return eval_result;
}

func pwd(nil)
/* DOCUMENT pwd
       -or- pwd()
     Prints out (subroutine form) or returns (function form) full path
     of current working directory.

   SEE ALSO: cd, lsdir. */
{
  if (! is_void(nil)) error, "unexpected non-nil argument";
  dir = cd(".");
  if (am_subroutine()) write, format="%s\n", dir;
  else return dir;
}

local unref; /* needed for documentation */
func __unref(&x) /* interpreted version */
/* DOCUMENT unref(x)
     returns X, destroying X in the process (useful to deal with temporary
     big arrays).  Written after Yorick's FAQ.
   SEE ALSO: eq_nocopy, swap. */
{
  local y;
  eq_nocopy, y, x;
  x = [];
  return y;
}
/* replace non-builtin function by interpreted one */
if (is_func(unref) != 2) unref =  __unref;

local swap; /* needed for documentation */
func __swap(&a, &b) /* interpreted version */
/* DOCUMENT swap, a, b;
     Exchanges  the contents  of variables  A and  B without  requiring any
     temporary copy.
   SEE ALSO: eq_nocopy, unref. */
{
  local tmp;
  eq_nocopy, tmp, a;
  eq_nocopy, a, b;
  eq_nocopy, b, tmp;
}
/* replace non-builtin function by interpreted one */
if (is_func(swap) != 2) swap =  __swap;

func grow_dimlist(&dimlist, arg)
/* DOCUMENT grow_dimlist, dimlist, next_argument
     builds a dimension  list DIMLIST, as used in  the array function.  Use
     like  this (all  extra arguments  in your  function are  considered as
     dimensions or dimension lists):
       func your_function(arg1, arg2, etc, ..)
       {
         dimlist = [0];
         while (more_args()) grow_dimlist, dimlist, next_arg();
         ...
       }
     After    this,   DIMLIST   will    be   an    array   of    the   form
     [ndims, dim1, dim2,...], compounded from the multiple arguments in the
     same way as the array function.  Another possibility is to define your
     function as:
       func your_function(arg1, arg2, etc, dimlist, ..)
       {
         while (more_args()) grow_dimlist, dimlist, next_arg();
         ...
       }
     But in this  latter case, if no DIMLIST  arguments given, DIMLIST will
     be [] instead of [0], which  will act the same in most situations.  If
     that possibility is unacceptible, you may add
       if (is_void(dimlist)) dimlist= [0];
     before/after the while loop.


   SEE ALSO: array, build_dimlist */
{
  if (is_array(arg)) {
    if (structof((n = arg(1)+0)) == long) {
      if (! (d1 = dimsof(arg)(1))) {
        if (is_void(dimlist)) {
          dimlist = [1, n];
        } else {
          grow, dimlist, n;
          ++dimlist(1);
        }
        return;
      } else if (d1 == 1) {
        if (is_void(dimlist)) {
          dimlist = long(arg);
          return;
        } else {
          if (n == numberof(arg)-1) {
            grow, dimlist, long(arg(2:0));
            dimlist(1) += n;
            return;
          }
        }
      }
      error, "bad dimension list";
    }
  } else if (is_void(arg)) {
    if (is_void(dimlist)) dimlist = [0];
    return;
  }
  error, "bad data type in dimension list";
}

func reform(x, ..)
/* DOCUMENT reform(x, dimlist, ...);
       -or- reform(x);
     returns array X reshaped according to dimension list DIMLIST.  Without
     DIMLIST, discards all dimensions equal to 1.  In most cases, prefer
     this to reshape.

   SEE ALSO: array, dimsof, grow_dimlist. */
{
  if (more_args()) {
    /* Build new dimension list. */
    dims = [0];
    while (more_args()) grow_dimlist, dims, next_arg();
  } else {
    /* Discard all dimensions equal to 1 */
    if (numberof(x) == 1) return x(1);
    dims = dimsof(x)(2:);
    dims = dims(where(dims>1));
    dims = [numberof(dims), dims];
  }
  y = array(structof(x), dims);
  y(*) = x(*);   /* will blow up if lengths differ */
  return y;
}

func undersample(a, nsub, which=, op=)
/* DOCUMENT undersample(a, nsub)
     Returns  array A  with all  (some)  dimensions divided  by NSUB.   The
     dimensions of interest must be a multiple of NSUB.

     Keyword  WHICH can  be used  to  specify the  dimension(s) to  shrink.
     Values  in  WHICH  less  or  equal  zero are  counted  from  the  last
     dimension.  By default, all dimensions of A get undersampled.

     Keyword OP can  be used to specify the range operator  to apply to the
     sets of NSUB adjacent values along the considered dimensions:
       OP=sum   to sum the values
       OP=avg   to average values
       OP=min   to keep the smallest value
       OP=max   to keep the largest value
     By default,  the median is  taken (WARNING: with the  median operator,
     the result depends in which order the dimensions of A are considered).

   SEE ALSO: median. */
{
  if (nsub < 1) error, "NSUB must be >= 1";
  if (nsub == 1) return a;
  if (! is_array(a)) error, "expecting an array";
  rank = dimsof(a)(1);
  if (is_void(which)) {
    which = indgen(rank);
  } else {
    which += (which <= 0)*rank;
    if (structof(which) != long)
      error, "bad data type for keyword WHICH";
    if (min(which) < 1 || max(which) > rank)
      error, "out of range dimension in keyword WHICH";
  }
  nw = numberof(which);
  noop = is_void(op); /* take median value */
  if (! noop && typeof(op) != "range") error, "OP must be nil or a range operator";
  dims = array(rank+1, rank+2);
  for (k=1 ; k<=nw ; ++k) {
    this = which(k);
    if (this != 1) a = transpose(a, [1,this]);
    dims(2:) = dimsof(a);
    dims(2) = nsub;
    dims(3) = dims(3)/nsub;
    (tmp = array(structof(a), dims))(*) = a(*);
    if (noop) {
      a = median(tmp, 1);
    } else {
      a = tmp(op,..);
    }
    tmp = []; /* free some memory */
    if (this != 1) a = transpose(a, [this,1]);
  }
  return a;
}

func spline_zoom(a, factor, rgb=)
/* DOCUMENT spline_zoom(a, fact)
     Return an array obtained by cubic spline interpolation of A with all
     its dimension multiplied by a factor FACT.  If keyword RGB is true the
     first dimsion of A is left unchanged.  If keyword RGB is not
     specified, then it is considered as true if A is a 3 dimensional array
     of 'char' with its first dimension equal to 3.
     
   SEE ALSO: spline, transpose. */
{
  if (! is_func(spline)) require, "spline.i";
  dims = dimsof(a);
  ndims = dims(1);

  if (is_void(rgb)) {
    /* RGB image? */
    rgb = (structof(a) == char && ndims == 3 && dims(2) == 3);
  }
  if (rgb) {
    a = transpose(a, 0);
    k = 1;
  } else {
    k = 0;
  }
  while (++k <= ndims) {
    dims = dimsof(a);
    n0 = dims(2);
    n1 = max(long(factor*n0 + 0.5), 1);
    if (n1 == 1) {
      a = a(avg,..)(..,-);
    } else {
      if (n1 != n0) {
        dims(2) = n1;
        b = array(double, dims);
        x0 = indgen(0:n0-1);
        x1 = indgen(0:n1-1)*((n0 - 1.0)/(n1 - 1.0));
        n = numberof(a)/n0;
        for (i=1 ; i<=n ; ++i) b(,i) = spline(a(,i), x0, x1);
        eq_nocopy, a, b;
      }
      if (ndims > 1) a = transpose(a, 0);
    }
  }
  if (rgb) return char(max(min(floor(a + 0.5), 255.0), 0.0));
  return a;
}

func map(__map__f, __map__x)
/* DOCUMENT map(f, input)
     Map scalar function F onto array/list argument INPUT to mimics
     element-wise unary operation.
   SEE ALSO: _map. */
{
  /* all locals here must have weird names, since the user's function may
     rely on external variables for arguments not varying in the source,
     or for accumulated outputs */
  if (is_array(__map__x)) {
    __map__y = array(structof((__map__i = __map__f(__map__x(1)))),
                     dimsof(__map__x));
    __map__y(1) = __map__i;
    __map__n = numberof(__map__x);
    for (__map__i=2 ; __map__i<=__map__n ; ++__map__i) {
      __map__y(__map__i) = __map__f(__map__x(__map__i));
    }
  } else if ((__map__n = typeof(__map__x)) == "list") {
    __map__y = __map__i = _lst(__map__f(_car(__map__x)));
    for (__map__x = _cdr(__map__x) ; ! is_void(__map__x) ;
         __map__x = _cdr(__map__x)) {
      _cat, __map__i, _lst(__map__f(_car(__map__x)));
      __map__i= _cdr(__map__i);
    }
  } else if (! is_void(__map__x)) {
    error, "unsupported data type \""+__map__n+"\"";
  }
  return __map__y;
}

/*---------------------------------------------------------------------------*/
/* STRING ROUTINES */

local _strlower, _strupper;
/* DOCUMENT local _strlower, _strupper;
     Private arrays to convert char to upper/lowercase letters.
   SEE ALSO strlower, strupper */
(_strlower = char(indgen(0:255)))(1+'A':1+'Z') = _strlower(1+'a':1+'z');
(_strupper = char(indgen(0:255)))(1+'a':1+'z') = _strupper(1+'A':1+'Z');

local strlower, strtolower; /* needed for documentation */
func __strlower(s) /* interpreted version */
/* DOCUMENT strlower(s)
       -or- strtolower(s)
     Convert a string or an array of strings S to lower case letters.
   SEE ALSO strupper */
{
  /* fool codger */ extern _strlower;
  n = numberof((r = array(string, dimsof(s))));
  for (i=1; i<=n; ++i) r(i)= string(&_strlower(1+*pointer(s(i))));
  return r;
}

local strupper, strtoupper; /* needed for documentation */
func __strupper(s) /* interpreted version */
/* DOCUMENT strupper(s)
       -or- strtoupper(s)
     Convert a string or an array of strings S to upper case letters.
   SEE ALSO strlower */
{
  /* fool codger */ extern _strupper;
  n = numberof((r = array(string, dimsof(s))));
  for (i=1; i<=n; ++i) r(i)= string(&_strupper(1+*pointer(s(i))));
  return r;
}

/* replace non-builtin functions by interpreted ones */
if (is_func(strupper) != 2) strupper = __strupper;
if (is_func(strlower) != 2) strlower = __strlower;
if (is_func(strtoupper) != 2) strtoupper = strupper;
if (is_func(strtoupper) != 2) strtolower = strtolower;

func strcut(str, len)
/* DOCUMENT strcut(str, len)
     Cut input scalar string STR in pieces of length less or equal LEN and
     return an array of such pieces.

   SEE ALSO strjoin */
{
  if ((str_len= strlen(str))<=len) return str;
  n= (str_len+len-1)/len;
  result= array(string, n);
  for (i=1, i1=1, i2=len ; i<=n ; ++i, i1+=len, i2+=len)
    result(i)= strpart(str, i1:i2);
  return result;
}

func strjoin(str, glue)
/* DOCUMENT strjoin(str)
       -or- strjoin(str, glue)
     Join strings from array STR into a single long string.  The string GLUE
     (default "") is used between each pair of element from STR.

   SEE ALSO strcut */
{
  if ((n= numberof(str)) >= 1) {
    s= str(1);
    if (glue) for (i=2 ; i<=n ; ++i) s+= glue+str(i);
    else      for (i=2 ; i<=n ; ++i) s+= str(i);
    return s;
  }
}

/*---------------------------------------------------------------------------*/
/* LOGICAL ROUTINES */
func _is_scalar(x) { return (is_array(x) && ! dimsof(x)(1)); }
/* DOCUMENT is_scalar(x)
     returns true if X is a scalar.

   SEE ALSO is_scalar, is_vector, is_matrix, is_array,
            is_integer, is_real, is_complex, is_numerical. */
if (!is_scalar) is_scalar = _is_scalar;

func _is_vector(x) { return (is_array(x) && dimsof(x)(1) == 1); }
/* DOCUMENT is_vector(x)
     returns true if array X is a vector (i.e. a 1D array).

   SEE ALSO is_scalar, is_vector, is_matrix, is_array,
            is_integer, is_real, is_complex, is_numerical. */
if (!is_vector) is_vector = _is_vector;

func _is_matrix(x) { return (is_array(x) && dimsof(x)(1) == 2); }
/* DOCUMENT is_matrix(x)
     returns true if array X is a matrix (i.e. a 2D array).

   SEE ALSO is_scalar, is_vector, is_matrix, is_array,
            is_integer, is_real, is_complex, is_numerical. */
if (!is_matrix) is_matrix = _is_matrix;


func _is_real(x)    { return ((s=structof(x))==double || s==float); }
/* DOCUMENT is_real(x)
     returns true if array X if of real type (i.e. double or float).

   SEE ALSO is_scalar, is_vector, is_matrix, is_array,
            is_integer, is_real, is_complex, is_numerical. */
if (!is_real) is_real = _is_real;


func _is_complex(x) { return structof(x)==complex; }
/* DOCUMENT is_integer(x)
     returns true if array X if of integer type.

   SEE ALSO is_scalar, is_vector, is_matrix, is_array,
            is_integer, is_real, is_complex, is_numerical. */
if (!is_complex) is_complex = _is_complex;


func _is_integer(x) 
/* DOCUMENT is_integer(x)
     returns true if array X if of integer type.

   SEE ALSO is_scalar, is_vector, is_matrix, is_array,
            is_integer, is_real, is_complex, is_numerical. */
{ return ((s=structof(x))==long || s==int || s==char || s==short);}
if (!is_integer) is_integer = _is_integer;


func _is_numerical(x) 
/* DOCUMENT is_numerical(x)
     returns true if array X if of numerical type.

   SEE ALSO is_scalar, is_vector, is_matrix, is_array,
            is_integer, is_real, is_complex, is_numerical. */
{ return ((s=structof(x))==long || s==double || s==int || s==char ||
          s==complex || s==short || s==float); }
if (!is_numerical) is_numerical = _is_numerical;


func _is_integer_scalar(x) 
/* DOCUMENT is_integer_scalar(x)
     Check whether or not X is an integer scalar.

   SEE ALSO is_scalar, is_integer. */
{ return (((s=structof(x))==long || s==int || s==short ||
           s==char) &&! dimsof(x)(1)); }
if (!is_integer_scalar) is_integer_scalar = _is_integer_scalar;

/*---------------------------------------------------------------------------*/
/* FILE ROUTINES */

func expand_file_name(_path)
/* DOCUMENT expand_file_name(path)
     Expand leading "~" in file name PATH which must be an array of
     strings (or a scalar string).
   SEE ALSO: strpart, cd, get_cwd, get_file_name, protect_file_name. */
{
  result = array(string, dimsof(_path));
  n = numberof(_path);
  cwd = get_cwd(); /* memorize current working directory */
  head = string(0);
  local path;
  for (i=1 ; i<=n ; ++i) {
    /* Expand leading "~" in file name(s). */
    eq_nocopy, path, _path(i);
    if (strpart(path, 1:1) == "~") {
      sread, path, format="%[^/]", head;
      home = cd(head);
      if (home) {
        cd, cwd; /* restore working directory */
        tail = strpart(path, strlen(head)+1:0);
        if (strlen(tail)) {
          /* remove trailing '/' in HOME to avoid two '/' in result */
          result(i) = strpart(home, 1:-1) + tail;
        } else if (strlen(home) > 1) {
          path = strpart(home, 1:-1);
        } else {
          eq_nocopy, path, home;
        }
      }
    }
    result(i) = path;
  }
  return result;
#if 0
  c = *pointer(path);
  i = numberof(c);
  expand = (i>1 && c(1) == '~');
  slash = (i>1 && c(-1) == '/');
  while (--i>0 && c(i) != '/')
    ;
#endif
}

func get_file_name(obj)
/* DOCUMENT get_file_name(obj)
     If OBJ is a stream, returns the path name of the file associated with
     the stream.  Otherwise if OBJ is an array of strings, expand leading
     "~" in the elements of OBJ.
     
   SEE ALSO: open, print, expand_file_name, protect_file_name. */
{
  /* Expand leading "~" if string object. */
  if (structof(obj) == string) {
    return expand_file_name(obj);
  }

  /* Check input and get description of stream by the print() command. */
  if ((id = typeof(obj)) == "stream") {
    id = 1;
    s = print(obj);
  } else if (id == "text_stream") {
    id = 2;
    s = print(obj)(2:);
  } else {
    error, "unexpected non-string, non-stream argument";
  }

  /* Join backslash terminated lines from print() result (another
     possibility would be to change the line length with `print_format' but
     there is no way to restore the previous line_lenght unles we building
     a wrapper around original `print_format' routine and make a
     substitution). */
  join = (strpart(s, 0:0) == "\\");
  if (anyof(join)) {
    r = array(string, (ns = numberof(s)) - sum(join) + join(0));
    i = j = 0;
    while (i < ns) {
      w = s(++i);
      while (join(i)) {
        w = strpart(w, :-1);
        if (++i > ns) break;
        w += s(i);
      }
      r(++j) = w;
    }
    s = r;
    w = r = [];
  }

  /* Recover the full path of the stream file from the joined lines. */
  if (id == 1) {
    /* Binary stream. */
    if (numberof(s) == 2) {
      w1 = w2 = string(0);
      if (sread(s(1), format="%[^:]", w1) == 1 &&
          sread(s(2), format="%[^/]", w2) == 1) {
        return strpart(s(2), strlen(w2)+1:0) + strpart(s(1), strlen(w1)+3:0);
      }
    }
    error, "unexpected binary stream descriptor";
  } else {
    /* Text stream. */
    if (numberof(s) == 1) {
      w = string(0);
      if (sread(s(1), format="%[^/]", w) == 1) {
        return strpart(s(1), strlen(w)+1:0);
      }
    }
    error, "unexpected text stream descriptor";
  }
}

local _protect_file_name_table, _protect_file_name_list;
func protect_file_name(path)
/* DOCUMENT protect_file_name(path)
     Protect special characters in PATH (apostrophes, quotes, $, *, ?, [,
     etc) by backslashes to avoid them being interpreted by the shell, for
     instance when using the system() builtin function.
   SEE ALSO: system, get_file_name, expand_file_name. */
{
  c = *pointer(path);
  n = numberof(c) - 1; /* same as strlen(path) */
  p = _protect_file_name_table(1 + c);
  r = array(char, 2*n + 1);
  i = 0;
  j = 0;
  while (++i <= n) {
    if (p(i)) r(++j) = '\\';
    r(++j) = c(i);
  }
  return string(&r);
}
_protect_file_name_list = ['$','&','!','#','?','*',
                           '[',']','{','}','<','>',
                           '\\','"','\'',
                           ' ','\t','\r','\n','\v','\f'];
(_protect_file_name_table = array(char, 256))(1 + _protect_file_name_list) = 1;

func read_ascii(file, compress=)
/* DOCUMENT read_ascii(file_or_name)
     Reads ascii numeric  data in columns from text  file.  FILE_OR_NAME is
     the name of the  file or an already open file stream.  The result is a
     NCOLUMNS-by-NLINES array of doubles.

     Data are  read as double values  arranged in columns  separated by any
     number of  spaces or tabs.  Comments  starting with a "#"  or an FILEy
     other character  which is not part of  a number are ignored  up to the
     end-of-line.  Blank lines  are ignored.  The first non-blank/commented
     line  gives the  number of  values per  column, for  subsequent lines.
     Subsequent lines  must have  the same number  of columns --  blanks in
     columns are  not permitted, use  0.0 instead.  However,  minimal error
     checking  is performed,  and if  the data  is not  really  in columns,
     read_ascii  can silently  fail to  interpret  your file  as you  would
     scanning it by eye.
     
     The  read operation will  be much  faster if  the number  of commented
     lines is  relatively small.   Blank lines cost  nothing, while  a line
     containing just a "#" is expensive.

     If the  file is specified by its  name, it may be  compressed in which
     case it get automatically decompressed while reading (see xopen).

   SEE ALSO: xopen, read. */
{
  /* open the file if it's not already open */
  if (structof(file) == string)
    file = xopen(file, compress=(is_void(compress) ? "auto" : compress));

  /* read lines one at a time until the "model" line which
   * determines the number of columns is discovered
   * assume the number of columns is less than 128 */
  x = array(0.0, 128);
  ncols = 0;
  while ((line = rdline(file))) {
    ncols = sread(line, x);
    if (ncols) break;          /* got a line with numbers */
  }
  if (! ncols) return [];

  nrows = 1;
  list = _lst([x(1:ncols)]);
  x = array(0.0, ncols, 10000/ncols + 1);
  for(;;) {
    /* try to grab at least 10000 numbers from the file
     * blank lines will be skipped, but any comments will
     * interrupt the read */
    if (! (n = read(file, x))) {
      /* if didn't get any, drop back to reading comments one
       * line at a time until we get some more numbers */
      while ((line = rdline(file))) {
	if ((n = sread(line, x))) break;
      }
      if (! line) break;    /* rdline detected end-of-file, n==0 too */
    }
    if (n%ncols) error, "data is not in columns";
    n /= ncols;

    /* grow the list the fast way, adding new values to its head
     * (adding to the tail would make growth an n^2 proposition,
     *  as would using the grow function) */
    list = _cat(x(,1:n), list);
    nrows += n;
  }

  /* pop chunks off list and reassemble result */
  x = array(0.0, ncols, nrows);
  for (i=nrows ; list ; list=_cdr(list)) {
    n = numberof(_car(list))/ncols;
    x(,i-n+1:i) = _car(list);
    i -= n;
  }

  return x;
}

func load_text(file, compress=)
/* DOCUMENT load_text(file)
     Returns all lines of text file as a vector of strings.  Returns nil if
     there are no lines to read.  FILE can be a file name or a text stream.
     In the latter case, the lines not yet read get returned.

     By  default,  if  the file  is  specified  by  its  name, it  will  be
     automatically decompressed if its  compressed; it is possible to force
     another  behaviour by  specifying a  value different  from  "auto" for
     keyword COMPRESS (see xopen).
     
   SEE ALSO: xopen, rdline, dump_text. */
{
  if (structof(file) == string)
    file = xopen(file, compress=(is_void(compress) ? "auto" : compress));
  text = rdline(file, 1000);
  while (text(0)) grow, text, rdline(file, numberof(text));
  if (is_array((j = where(text)))) return text(j);
}

func dump_text(file, text, compress=, level=, preserve=)
/* DOCUMENT dump_text, file, text;
     Dump every  elements of string array  TEXT as individual  lines into a
     text file.  FILE can be a file name or a text stream.

     If the file is specified by  its name, keywords COMPRESS, LEVEL can be
     used to specify  a compression method (see xopen).   The default value
     of COMPRESS is  "auto" (i.e., method guessed on the  basis of the file
     extension).

     If the file  is specified by its name, keyword PRESERVE  can be set to
     true to avoid overwriting an existing file.
     
   SEE ALSO: xopen, rdline, load_text. */
  
{
  if ((anything = ! is_void(text)) && structof(text) != string)
    error, "expecting array of strings or nil";
  if (structof(file) == string) {
    file = xopen(file, "w", compress=(is_void(compress) ? "auto" : compress),
                 level=level, preserve=preserve);
  }
  if (anything) write, file, format="%s\n", text;
}

func guess_compression(filename)
/* DOCUMENT guess_compression, filename;
       -or- guess_compression(filename)  
     Guess  which compression  program was  used to  produce  file FILENAME
     according to first  bytes of this file.  When  called as a subroutine,
     the file name is printed out  with the name of the compression program
     if any.  If called as a function, the result is an integer:
       1 - if file compressed with "gzip";
       2 - if file compressed with "bzip2";
       3 - if file compressed with "compress";
       4 - if file compressed with "pack";
       0 - otherwise.

  SEE ALSO: xopen. */
{
  /* according to information in /etc/magic:
   *
   * method    min.            bytes
   *           size      0    1    2    3
   * --------- ----    ---- ---- ---- ----
   * pack:       3     \037 \036
   * compress:   3     \037 \235
   * gzip:      20     \037 \213   c         (1)
   * bzip2:     14      'B'  'Z'  'h'   c    (2)
   *
   *   (1)  if c<8, compression level, else if c==8 deflated
   *   (2)  with '0' <= c <= '9', block size = c*100kB
   *   minimum size has been computed for an empty file:
   *     pack     ->  3 bytes
   *     compress ->  3 bytes
   *     bzip2    -> 14 bytes
   *     gzip     -> 20 bytes, if compression level is specified
   *                 24 bytes, otherwise
   */
  magic = array(char, 4);
  n = _read(open(filename, "rb"), 0, magic);
  if ((c = magic(1)) == '\037') {
    if ((c = magic(2)) == '\213') {
      if (! am_subroutine()) return 1; /* gzip */
      write, format="%s: compressed with \"gzip\"\n", filename;
      return;
    } else if (c == '\235') {
      if (! am_subroutine()) return 3; /* compress */
      write, format="%s: compressed with \"compress\"\n", filename;
      return;
    } else if (c == '\036') {
      if (! am_subroutine()) return 4; /* pack */
      write, format="%s: compressed with \"pack\"\n", filename;
      return;
    }   
  } else if (c == 'B' && magic(2) == 'Z' && magic(3) == 'h' &&
             '0' <= (c = magic(4)) && c <= '9') {
    if (! am_subroutine()) return 2; /* bzip2 */
    write, format="%s: compressed with \"bzip2\" (block size = %d kB)\n",
      filename, 100*(c - '0');
    return;
  }
  if (! am_subroutine()) return 0;
  write, format="%s: uncompressed?\n", filename;
}

func xopen(filename, filemode, preserve=, nolog=, compress=, level=, prims=)
/* DOCUMENT xopen(filename)
       -or- xopen(filename, filemode)
     Opens the file FILENAME according to FILEMODE (both are strings).  The
     return value is an IOStream (or just stream for short).  When the last
     reference to this return value is discarded, the file will be closed.
     The file can also be  explicitly closed with the close function (which
     see).  The  FILEMODE (default  "r" -- open  an existing text  file for
     reading) determines whether  the file is to be  opened in read, write,
     or update mode,  and whether writes are restricted  to the end-of-file
     (append mode).  FILEMODE also determines whether the file is opened as
     a text  file or  as a  binary file.  FILEMODE  can have  the following
     values, which are the same as for the ANSI standard fopen function:
         "r"   - read only
         "w"   - write only, random access, existing file overwritten
         "a"   - write only, forced to end-of-file, existing file preserved
         "r+"  - read/write, random access, existing file preserved
         "w+"  - read/write, random access, existing file overwritten
         "a+"  - read/write, reads random access, writes forced to
                 end-of-file, existing file preserved
         "rb"  "wb"  "ab"  "r+b"  "rb+"  "w+b"  "wb+"  "a+b"  "ab+"
                 without b means text file, with b means binary file
     
     Keyword COMPRESS can be used  to specify compression method for a text
     file open for reading (FILEMODE="r") or writing (FILEMODE="w") only --
     (de)compression  is unsupported in  append mode  or for  binary files.
     The value of keyword COMPRESS can be a scalar string or an integer:
          "auto"     - guess compression according to first bytes of file
                       in read  mode, or according to file extension in 
                       write mode: ".gz" for gzip, ".bz2" for bzip2 and
                       ".Z" for compress.
       0  "none"     - no (de)compression
       1  "gzip"     - use gzip to (de)compress
       2  "bzip2"    - use bzip2 to (de)compress
       3  "compress" - use compress to (de)compress
       4  "pack"     - use pack to (de)compress
     The default  value for COMPRESS is  "auto" in read mode  and "none" in
     write mode.  Note that "gzip", "bzip2", "pack" and "compress" commands
     must  exists in  your  $PATH for  compressing  with the  corresponding
     methods.  Decompression of files compressed with "pack" and "compress"
     is  done  by  "gzip".   If  keyword  COMPRESS  is  explicitely  0,  no
     decompression  is ever  applied;  if keyword  COMPRESS is  explicitely
     non-zero, the  file must have been compressed.   The compression level
     for gzip  and bzip2  can be  specified as an  integer value  thanks to
     keyword LEVEL.
       
     Keyword PRIMS  can be used  to specify primitives data  type different
     than  the native  ones for  binary files  (PRIMS is  ignored  for text
     files).  PRIMS can  be a scalar string (i.e.,  "alpha", "mac", "sun3",
     "cray", "macl", "sun", "dec",  "pc", "vax", "vaxg", "i86", "sgi64", or
     "xdr") or  a 32-element  vector of long's  as taken  by set_primitives
     (which see).

     When a binary file is created, it is possible to avoid the creation of
     the log-file FILENAME+"L" by setting keyword NOLOG to true.

     Keyword PRESERVE can  be set to true to  avoid overwriting an existing
     file when FILENAME is open for writing (i.e. with a "w" in FILEMODE).

     
   BUGS:
     If (de)compression is used, FILENAME must not contain any double quote
     character (").

     
   SEE ALSO: close, guess_compression, open, popen, set_primitives. */
{
  if (is_void(filemode) || filemode == "r") {
    /* Open file for reading in text mode. */
    compress = __xopen_get_compress(compress, filename, 1);
    if (! compress) return open(filename, "r");
    if (compress == 2) return popen("bzip2 -dc \"" + filename + "\"", 0);
    if (compress == -1) error, "bad value for keyword COMPRESS";
    return popen("gzip -dc \"" + filename + "\"", 0);
  }

  if (preserve && strmatch(filemode, "w") && open(filename, "r", 1))
    error, "file \""+filename+"\" already exists";

  filemode
  
  if (filemode == "w") {
    /* Open file for writing in text mode. */
    compress = __xopen_get_compress(compress, filename, 0);
    if (! compress) return open(filename, filemode);
    if (compress == 1) {
      if (is_void(level)) command = swrite(format="gzip > \"%s\"", filename);
      else command = swrite(format="gzip -%d > \"%s\"", level, filename);
    } else if (compress == 2) {
      if (is_void(level)) command = swrite(format="bzip2 > \"%s\"", filename);
      else command = swrite(format="bzip2 -%d > \"%s\"", level, filename);
    } else if (compress == 3) {
      command = swrite(format="compress > \"%s\"", filename);
    } else if (compress == 2) {
      command = swrite(format="pack > \"%s\"", filename);
    } else {
      error, "bad value for keyword COMPRESS";
    }
    command
    return popen(command, 1);
  }

  /* Open file for other modes. */
  if (! (is_void(compress) || compress == 0 || compress == "none"))
    error, "(de)compression unsupported in mode \""+filemode+"\"";
  if (binary && nolog) {
    /* will remove log-file unless it already exists */
    logfile = filename + "L";
    if (open(logfile, "r", 1)) logfile = [];
  } else {
    logfile = [];
  }
  file = open(filename, filemode);
  if (logfile) remove, logfile;

  /* Return open file after optionally installing primitive data types. */
  if (is_void(prims) || ! strmatch(filemode, "b")) return file;
  if ((s = structof(prims)) == string) {
    if (! dimsof(prims)(1)) {
      if (prims != "set" && prims != "get" &&
          is_func((sym = symbol_def(prims+"_primitives"))) == 1) {
        sym, file;
        return file;
      } else if (is_array((sym = symbol_def("__"+prims))) &&
                 structof(sym) == long && numberof(sym) == 32 &&
                 dimsof(sym)(1) == 1) {
        set_primitives, file, sym;
        return file;
      }
    } else if (s == long && numberof(prims) == 32 && dimsof(prims)(1) == 1) {
      set_primitives, file, prims;
      return file;
    }
  }
  error, "bad value for keyword PRIMS";
}

func __xopen_get_compress(compress, filename, for_reading)
/* DOCUMENT __xopen_get_compress(compress, filename, for_reading)
     Private function called by xopen.

   SEE ALSO xopen. */
{
  if (is_void(compress)) {
    if (for_reading) return guess_compression(filename);
    return 0;
  } else if (is_array(compress) && ! dimsof(compress)(1)) {
    if ((s = structof(compress)) == string) {
      if (compress == "auto") {
        if (for_reading) return guess_compression(filename);
        if (strpart(filename, -2:0) == ".gz" ) return 1; /* gzip */
        if (strpart(filename, -3:0) == ".bz2") return 2; /* bzip2 */
        if (strpart(filename, -1:0) == ".Z"  ) return 3; /* compress */
        return 0;
      }
      if (compress == "none"    ) return 0;
      if (compress == "gzip"    ) return 1;
      if (compress == "bzip2"   ) return 2;
      if (compress == "compress") return 3;
      if (compress == "pack"    ) return 4;
    } else if ((s == long || s == int || s == short || s == char) &&
               compress >= 0 && compress <= 4) {
      return compress;
    }
  }
  return -1;
}


/*---------------------------------------------------------------------------*/
/* PDB FILES */

func pdb_list(file)
/* DOCUMENT pdb_list, file;
       -or- pdb_list(file)
     Lists contents of PDB binary file.  FILE can be either a file name or
     a binary stream.
     
   SEE ALSO: createb, openb, restore, pdb_restore_all. */
{
  if (structof(file) == string) file = openb(file);
  vars = get_vars(file);
  if (! am_subroutine()) return vars;
  title = ["Non-record variables", "    Record variables"];
  for (i=1 ; i<=2 ; ++i) {
    write, format="%s:", title(i);
    if (numberof(*vars(i))) {
      write, format=" %s", *vars(i);
      write, format="%s\n", "";
    } else {
      write, format="%s\n", " <none>";
    }
  }
}

func pdb_restore_all(_f_i_l_e_)
/* DOCUMENT pdb_restore_all, file;
     Restore all non-record variables of a PDB file.  FILE can be either a
     file name or a binary stream. 
     
   SEE ALSO: createb, openb, restore, pdb_list. */
{
  if (structof(_f_i_l_e_) == string) _f_i_l_e_ = openb(_f_i_l_e_);
  restore, _f_i_l_e_;
}

/*---------------------------------------------------------------------------*/
/* PROFILING ROUTINES */

local _timer_stamp;
func timer_start {
  extern _timer_stamp;
  _timer_stamp = array(double, 3);
  timer, _timer_stamp;
}
func timer_elapsed(count)
/* DOCUMENT timer_start;
       -or- timer_elapsed;
       -or- timer_elapsed, count;
       -or- timer_elapsed()
       -or- timer_elapsed(count)
     The  subroutine  timer_start (re)starts  the  timer  and the  function
     timer_elapsed   computes  the   elapsed  time   since  last   call  to
     timer_start.  If COUNT is given, the elapsed time is divided by COUNT.
     When  called as  a subroutine,  timer_elapsed prints  out  the elapsed
     time; when  called as a  function, it returns  [CPU,SYSTEM,WALL], with
     all three  times measured in seconds.   The two functions  make use of
     external variable _timer_stamp to memorize the initiale times.

     For instance:
       timer_start;
       ...             // some code to be profiled
       timer_elapsed;
       ...             // some more code to be profiled
       timer_elapsed;  // prints out _total_ elapsed time

  SEE ALSO: timer. */
{
  extern _timer_stamp;
  elapsed = _timer_stamp;
  timer, elapsed;
  elapsed -= _timer_stamp;
  if (! is_void(count)) elapsed /= count;
  if (am_subroutine()) {
    write, format="cpu=%g, system=%g, wall=%g\n",
      elapsed(1), elapsed(2), elapsed(3);
  } else {
    return elapsed;
  }
}

/*---------------------------------------------------------------------------*/

func _stat_worker(x)
/* DOCUMENT _stat_worker(x)
     Private routine used by stat, returns vector of double's:
        [min(X), max(X), avg(X), std(X)]
     where std(X) is the standard deviation of X.

   SEE ALSO stat. */
{
  if (structof(x)!=double) x= double(x);
  avg_x= avg(x);
  dx = x - avg_x;
  return [min(x), max(x), avg_x, sqrt(avg(dx*dx))];
}

func stat(..)
/* DOCUMENT stat, x, ...
     Print out statistics and information for all the arguments. */
{
  ith= 0;
  while (more_args()) {
    ++ith;
    x= next_arg();
    write, format="%2d: ", ith;
    if (is_array(x)) {
      write, format="array(%s", typeof(x);
      dims= dimsof(x);
      n= numberof(dims);
      for (k=2 ; k<=n ; ++k) write, format=",%d", dims(k);
      type= structof(x);
      is_numerical= (type==double || type==long || type==int || type==char ||
                     type==complex || type==float || type==short);
      write, format=")%s", (is_numerical ? " " : "\n");
      if (is_numerical) {
        fmt= "min=%g max=%g avg=%g std=%g\n";
        if (type==complex) {
          s= _stat_worker(double(x));
          write, format="\n         real part: "+fmt, s(1), s(2), s(3), s(4);
          s= _stat_worker(x.im);
          write, format="    imaginary part: "+fmt, s(1), s(2), s(3), s(4);
          s= _stat_worker(abs(x));
          write, format="           modulus: "+fmt, s(1), s(2), s(3), s(4);
        } else {
          s= _stat_worker(x);
          write, format=fmt, s(1), s(2), s(3), s(4);
        }
      }
    } else {
      write, format="%s, %s\n", typeof(x), strjoin(print(x));
    }
  }
}

/*---------------------------------------------------------------------------*/

func smooth(a, level)
/* DOCUMENT smooth(a)
       -or- smooth(a, level)
     Returns array A smoothed along its dimensions.  I.e. for a 1D array:
       smooth(A) = A(pcen)(zcen)
     for a 2D array:
       smooth(A) = A(pcen,pcen)(zcen,zcen)
     ... (up to 6 dimensions).

     For a greater number of dimensions,  each  direction  is  smoothed and
     transposed in turn: apart from rounding errors, the result is the same
     but the computation time is approximately  multiplied  by  3.   If you
     oftenly smooth arrays with more than 6 dimensions you may  think about
     modifying the source...

     Optional argument  LEVEL  (default  1)  set  the  number  of  time the
     smoothing operation is performed.

   PROPERTIES OF THE SMOOTHING OPERATOR:
     (i)   The smoothing operator is linear and symmetric.  For instance,
           for a vector, A, smooth(A)=S(,+)*A(+) where the matrix S is
           tridiagonal:
                    [3 1         ]
                    [1 2 1       ]
                    [  1 2 1     ]
             0.25 * [   \ \ \    ]    where, to improve readability,
                    [    \ \ \   ]    missing values are all zero.
                    [     1 2 1  ]
                    [       1 2 1]
                    [         1 3]
           You can, in principle, reverse the smoothing operation with
           TDsolve along each dimensions of smooth(A).  Note:For a vector
           A, the operator S-I applied to A (where I is the identity
           matrix) is the finite difference 2nd derivatives of A (but for
           the edges).

     (ii)  The smoothing operator does not change the sum of the element
           values of its argument, i.e.: sum(smooth(A)) = sum(A).

     (iii) Only an array with all elements having the same value is
           invariant by the smoothing operator.  In fact "slopes" along
           dimensions of A are almost invariant, only the values along the
           edges are changed.

     The symmetry of the smoothing operator is important for the
     computation of gradients.  For instance, let Y = smooth(X) and DQ_DY
     be the gradient of a scalar function Q with respect to Y, then the
     gradient of Q with respect to X is simply: DQ_DX = smooth(DQ_DY)

   TO DO:
     By default A is smoothed along all its dimensions, but the list
     of dimensions to smooth can be specified with keyword WHICH.  As
     usual, negative dimensions are taken as offset from the last one.

     If keyword WRAP is true (non-nil and non-zero) a wrapped version
     of the operator (with same properties but no longer tridiagonal)
     is applied instead.  This is suitable for periodic arrays (e.g.
     FFT transformed arrays).

   SEE ALSO: TDsolve. */
{
  n= dimsof(a)(1);
  if (is_void(level) || level == 1) {
    if (n == 1)
      return a(pcen)(zcen);
    if (n == 2)
      return a(pcen,pcen)(zcen,zcen);
    if (n == 3)
      return a(pcen,pcen,pcen)(zcen,zcen,zcen);
    if (n == 4)
      return a(pcen,pcen,pcen,pcen)(zcen,zcen,zcen,zcen);
    if (n == 5)
      return a(pcen,pcen,pcen,pcen,pcen)(zcen,zcen,zcen,zcen,zcen);
    if (n == 6)
      return a(pcen,pcen,pcen,pcen,pcen,pcen)(zcen,zcen,zcen,zcen,zcen,zcen);
    while (n--)
      a= transpose(a(pcen,..)(zcen,..));
    return a;
  }
  if (n == 1) {
    for (i=1; i<=level; i++)
      a= a(pcen)(zcen);
  } else if (n == 2) {
    for (i=1; i<=level; i++)
      a= a(pcen,pcen)(zcen,zcen);
  } else if (n == 3) {
    for (i=1; i<=level; i++)
      a= a(pcen,pcen,pcen)(zcen,zcen,zcen);
  } else if (n == 4) {
    for (i=1; i<=level; i++)
      a= a(pcen,pcen,pcen,pcen)(zcen,zcen,zcen,zcen);
  } else if (n == 5) {
    for (i=1; i<=level; i++)
      a= a(pcen,pcen,pcen,pcen,pcen)(zcen,zcen,zcen,zcen,zcen);
  } else if (n == 6) {
    for (i=1; i<=level; i++)
      a= a(pcen,pcen,pcen,pcen,pcen,pcen)(zcen,zcen,zcen,zcen,zcen,zcen);
  } else {
    while (n--) {
      for (i=1; i<=level; i++)
	a= a(pcen,..)(zcen,..);
      a= transpose(a);
    }
  }
  return a;
}

/*---------------------------------------------------------------------------*/