n行ごとの平均をとる

スペクトル的なデータを取ったものの各ビンのデータ数が小さくて統計が足りないのでいくつかの行を平均してやってまとめたいというときはgnuplotのsmooth uniqをつかったりしていた*1のだけど、これは特にlog scaleでは大いにbuggy*2なのでスクリプトを作ってみた。需要あるとは思えないがclispスクリプトを書く練習として。

#channel data1	data2
1	2	2000
2	4	4000
3	4	8000
4	5	9000
5	6	10000
6	8	10000
(snip)
96	1	1000
97	0	500
98	1	200
99	0	100
100	0	100

みたいなデータ data.dat があった時に

$ clisp 2 data.dat

とすると data.dat.avr として 2行ごとの平均を各列についてとったものを出力する:

## averaged every 2 lines
#channel data1	data2
1.5 3.0 3000.0 
3.5 4.5 8500.0 
5.5 7.0 10000.0 

96.5 0.5 750.0 
98.5 0.5 150.0 
100.0 0.0 100.0 

というもの。空行とコメント行(#で始まる行)はそのまま出力する。ファイルは複数指定可。

#!/usr/bin/clisp
;; usage
;; ./average-lines COUNT file1 file2 ... fileN

;; averages COUNT lines together, each columns
;; for gnuplot data file

(defun my-command-line ()
  (or 
   #+SBCL *posix-argv*  
   #+LISPWORKS system:*line-arguments-list*
   #+CMU extensions:*command-line-words*
   #+CLISP *args*
   #+Allegro (sys:command-line-arguments)
   nil)) ;http://d.hatena.ne.jp/cl-intern/20070809#p2

(defmacro while (pred &body body)
  `(loop (progn (unless ,pred (return)) ,@body)))

(defvar n (read-from-string (car (my-command-line))))
(unless (typep n 'integer) (error "1st arg must be positive integer"))

(dolist (f (cdr (my-command-line)))
  (let ((count 0) (sums ()) buf)
    (with-open-file (is f :direction :input)
      (with-open-file (os (concatenate 'string f ".avr") :direction :output :if-exists :rename-and-delete :if-does-not-exist :create)
	(format os "## averaged every ~D lines~%" n)
	(labels ((output-avrs (sums) (format os "~{~F ~}~%" (mapcar #'(lambda (x) (/ x count)) sums))))
	  (while (setf buf (read-line is nil nil))
	    (cond
	      ((or (zerop (length buf)) (char= (aref buf 0) #\Return)) ;blank line: separetes data blocks
	       (when (not (zerop count))
		 (output-avrs sums))
	       (setf sums ()) (setf count 0)
	       (terpri os))
	      ((char= #\# (aref buf 0))	; comment line: output the comment line without change
	       (write-line buf os))
	      (t			; ordinary data line
	       (if (null sums)
		   (setf sums (let ((lst))
				(loop
				   (multiple-value-bind (val pos) (read-from-string buf nil (values nil 0))
				     (setf buf (subseq buf pos))
				     (if val (push val lst) (return))))
				(nreverse lst)))
		   (dotimes (i (length sums))
		     (multiple-value-bind (val pos) (read-from-string buf)
		       (incf (nth i sums) val)
		       (setf buf (subseq buf pos)))))
	       (incf count)
	       (when (>= count n)
		 (output-avrs sums)
		 (setf sums ()) (setf count 0)))))
	  (when (not (zerop count))
	    (output-avrs sums)))))))