[lisplab-cvs] r197 - in trunk: . src/io
Jørn Inge Vestgården
jivestgarden at common-lisp.net
Sat Dec 4 10:04:01 UTC 2010
Author: jivestgarden
Date: Sat Dec 4 05:04:00 2010
New Revision: 197
Log:
Added simple save and load utility
Added:
trunk/src/io/ieee-floats.lisp
trunk/src/io/saveload.lisp
Modified:
trunk/lisplab.asd
trunk/package.lisp
Modified: trunk/lisplab.asd
==============================================================================
--- trunk/lisplab.asd (original)
+++ trunk/lisplab.asd Sat Dec 4 05:04:00 2010
@@ -130,11 +130,14 @@
;;
;; IO (level 3)
;;
- (:module :src/io
+ (:module :src/io
:depends-on (:src/matrix2)
+ :serial t
:components
((:file "level3-io-interface")
- (:file "level3-io")))
+ (:file "level3-io")
+ (:file "ieee-floats")
+ (:file "saveload")))
;;
;; Linear algebra lisp implementation (Level 3)
Modified: trunk/package.lisp
==============================================================================
--- trunk/package.lisp (original)
+++ trunk/package.lisp Sat Dec 4 05:04:00 2010
@@ -55,6 +55,12 @@
"INIT-THREADS"
"CLEANUP-THREADS"
+ ;; Infix notation
+ "*SEPARATORS*"
+ "W/INFIX"
+ "INFIX->PREFIX"
+ "PREFIX->INFIX"
+
;; Numerical constants
"%I"
"%E"
@@ -124,12 +130,6 @@
".ERFC"
".GAMMA"
- ;; Infix notation
- "*SEPARATORS*"
- "W/INFIX"
- "INFIX->PREFIX"
- "PREFIX->INFIX"
-
;; Now the matrix stuff
;; Matrix classes
"MATRIX-BASE"
@@ -214,6 +214,8 @@
"PSWRITE"
"DLMREAD"
"DLMWRITE"
+ "MSAVE"
+ "MLOAD"
;; ODE solvers
"EULER" ; todo change name
Added: trunk/src/io/ieee-floats.lisp
==============================================================================
--- (empty file)
+++ trunk/src/io/ieee-floats.lisp Sat Dec 4 05:04:00 2010
@@ -0,0 +1,137 @@
+;;; Functions for converting floating point numbers represented in
+;;; IEEE 754 style to lisp numbers.
+;;;
+;;; See http://common-lisp.net/project/ieee-floats/
+
+(defpackage :ieee-floats
+ (:use :common-lisp)
+ (:export :make-float-converters
+ :encode-float32
+ :decode-float32
+ :encode-float64
+ :decode-float64))
+
+(in-package :ieee-floats)
+
+;; The following macro may look a bit overcomplicated to the casual
+;; reader. The main culprit is the fact that NaN and infinity can be
+;; optionally included, which adds a bunch of conditional parts.
+;;
+;; Assuming you already know more or less how floating point numbers
+;; are typically represented, I'll try to elaborate a bit on the more
+;; confusing parts, as marked by letters:
+;;
+;; (A) Exponents in IEEE floats are offset by half their range, for
+;; example with 8 exponent bits a number with exponent 2 has 129
+;; stored in its exponent field.
+;;
+;; (B) The maximum possible exponent is reserved for special cases
+;; (NaN, infinity).
+;;
+;; (C) If the exponent fits in the exponent-bits, we have to adjust
+;; the significand for the hidden bit. Because decode-float will
+;; return a significand between 0 and 1, and we want one between 1
+;; and 2 to be able to hide the hidden bit, we double it and then
+;; subtract one (the hidden bit) before converting it to integer
+;; representation (to adjust for this, 1 is subtracted from the
+;; exponent earlier). When the exponent is too small, we set it to
+;; zero (meaning no hidden bit, exponent of 1), and adjust the
+;; significand downward to compensate for this.
+;;
+;; (D) Here the hidden bit is added. When the exponent is 0, there is
+;; no hidden bit, and the exponent is interpreted as 1.
+;;
+;; (E) Here the exponent offset is subtracted, but also an extra
+;; factor to account for the fact that the bits stored in the
+;; significand are supposed to come after the 'decimal dot'.
+
+(defmacro make-float-converters (encoder-name
+ decoder-name
+ exponent-bits
+ significand-bits
+ support-nan-and-infinity-p)
+ "Writes an encoder and decoder function for floating point
+numbers with the given amount of exponent and significand
+bits (plus an extra sign bit). If support-nan-and-infinity-p is
+true, the decoders will also understand these special cases. NaN
+is represented as :not-a-number, and the infinities as
+:positive-infinity and :negative-infinity. Note that this means
+that the in- or output of these functions is not just floating
+point numbers anymore, but also keywords."
+ (let* ((total-bits (+ 1 exponent-bits significand-bits))
+ (exponent-offset (1- (expt 2 (1- exponent-bits)))) ; (A)
+ (sign-part `(ldb (byte 1 ,(1- total-bits)) bits))
+ (exponent-part `(ldb (byte ,exponent-bits ,significand-bits) bits))
+ (significand-part `(ldb (byte ,significand-bits 0) bits))
+ (nan support-nan-and-infinity-p)
+ (max-exponent (1- (expt 2 exponent-bits)))) ; (B)
+ `(progn
+ (defun ,encoder-name (float)
+ ,@(unless nan `((declare (type float float))))
+ (multiple-value-bind (sign significand exponent)
+ (cond ,@(when nan `(((eq float :not-a-number)
+ (values 0 1 ,max-exponent))
+ ((eq float :positive-infinity)
+ (values 0 0 ,max-exponent))
+ ((eq float :negative-infinity)
+ (values 1 0 ,max-exponent))))
+ ((zerop float)
+ (values 0 0 0))
+ (t
+ (multiple-value-bind (significand exponent sign) (decode-float float)
+ (let ((exponent (+ (1- exponent) ,exponent-offset))
+ (sign (if (= sign 1.0) 0 1)))
+ (unless (< exponent ,(expt 2 exponent-bits))
+ (error "Floating point overflow when encoding ~A." float))
+ (if (< exponent 0) ; (C)
+ (values sign (ash (round (* ,(expt 2 significand-bits) significand)) exponent) 0)
+ (values sign (round (* ,(expt 2 significand-bits) (1- (* significand 2)))) exponent))))))
+ (let ((bits 0))
+ (declare (type (unsigned-byte ,total-bits) bits))
+ (setf ,sign-part sign
+ ,exponent-part exponent
+ ,significand-part significand)
+ bits)))
+
+ (defun ,decoder-name (bits)
+ (declare (type (unsigned-byte ,total-bits) bits))
+ (let* ((sign ,sign-part)
+ (exponent ,exponent-part)
+ (significand ,significand-part))
+ ,@(when nan `((when (= exponent ,max-exponent)
+ (return-from ,decoder-name
+ (cond ((not (zerop significand)) :not-a-number)
+ ((zerop sign) :positive-infinity)
+ (t :negative-infinity))))))
+ (if (zerop exponent) ; (D)
+ (setf exponent 1)
+ (setf (ldb (byte 1 ,significand-bits) significand) 1))
+ (unless (zerop sign)
+ (setf significand (- significand)))
+ (scale-float (float significand ,(if (> total-bits 32) 1.0d0 1.0))
+ (- exponent ,(+ exponent-offset significand-bits)))))))) ; (E)
+
+;; And instances of the above for the common forms of floats.
+(make-float-converters encode-float32 decode-float32 8 23 nil)
+(make-float-converters encode-float64 decode-float64 11 52 nil)
+
+;;; Copyright (c) 2006 Marijn Haverbeke
+;;;
+;;; This software is provided 'as-is', without any express or implied
+;;; warranty. In no event will the authors be held liable for any
+;;; damages arising from the use of this software.
+;;;
+;;; Permission is granted to anyone to use this software for any
+;;; purpose, including commercial applications, and to alter it and
+;;; redistribute it freely, subject to the following restrictions:
+;;;
+;;; 1. The origin of this software must not be misrepresented; you must
+;;; not claim that you wrote the original software. If you use this
+;;; software in a product, an acknowledgment in the product
+;;; documentation would be appreciated but is not required.
+;;;
+;;; 2. Altered source versions must be plainly marked as such, and must
+;;; not be misrepresented as being the original software.
+;;;
+;;; 3. This notice may not be removed or altered from any source
+;;; distribution.
Added: trunk/src/io/saveload.lisp
==============================================================================
--- (empty file)
+++ trunk/src/io/saveload.lisp Sat Dec 4 05:04:00 2010
@@ -0,0 +1,160 @@
+;;; Lisplab, saveload.lisp
+;;; Input output operations in lisplab-specific format
+
+;;; Copyright (C) 2009 Joern Inge Vestgaarden
+;;;
+;;; 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.
+;;;
+;;; You should have received a copy of the GNU General Public License along
+;;; with this program; if not, write to the Free Software Foundation, Inc.,
+;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+;;; Some simple home made protocol for saving and loading matrices.
+;;; So far it only works for double float matrices.
+;;;
+;;; The file format is simply
+;;;
+;;; (nonce type header-length (metadata) rows cols (data))
+;;;
+;;; nonce, type, header-length, rows, and cols are 32 bit unsigned integers
+;;; the rest of the hearder is for metadata and is currently skipped
+;;; The data is 64 bit floats, in row major order
+;;; (and hopefully the numbers are stored as ieee compatible, big endian.)
+
+;;; In principle, lisplab should store and save matrices in some standard data format,
+;;; but thats a lot of work to implement.
+
+
+(in-package :lisplab)
+
+;;;; First some general stuff
+
+(defgeneric msave (stream-or-file matrix)
+ (:documentation "Writes the matrix in a lisplab-specific format."))
+(defgeneric mload (stream-or-file)
+ (:documentation "Reads a matrix coded in a lisplab-specific format."))
+
+(defmethod msave ((name pathname)
+ (a matrix-base))
+ (with-open-file (stream name
+ :direction :output
+ :if-exists :supersede
+ :element-type 'unsigned-byte)
+ (msave stream a)))
+
+(defmethod msave ((name string)
+ (a matrix-base))
+ (msave (pathname name) a))
+
+(defmethod mload ((name pathname))
+ (with-open-file (stream name
+ :direction :input
+ :element-type 'unsigned-byte)
+ (mload stream)))
+
+(defmethod mload ((name string))
+ (mload (pathname name)))
+
+;;;; Some helper functions
+
+(defun encode-ui32 (a i off)
+ "Writes four bytes to the byte array in big endian format."
+ (setf (aref a off) (ldb (byte 8 24) i)
+ (aref a (+ off 1)) (ldb (byte 8 16) i)
+ (aref a (+ off 2)) (ldb (byte 8 8) i)
+ (aref a (+ off 3)) (ldb (byte 8 0) i)))
+
+(defun decode-ui32 (a off)
+ "Reads a four byte integer from the byte array in big endian format."
+ (logior (ash (aref a off) 24)
+ (ash (aref a (+ off 1)) 16)
+ (ash (aref a (+ off 2)) 8)
+ (aref a (+ off 3))))
+
+(defun encode-ui64 (a i off)
+ "Writes eight bytes to the byte array in big endian format."
+ (setf (aref a off) (ldb '(8 . 56) i)
+ (aref a (+ off 1)) (ldb '(8 . 48) i)
+ (aref a (+ off 2)) (ldb '(8 . 40) i)
+ (aref a (+ off 3)) (ldb '(8 . 32) i)
+ (aref a (+ off 4)) (ldb '(8 . 24) i)
+ (aref a (+ off 5)) (ldb '(8 . 16) i)
+ (aref a (+ off 6)) (ldb '(8 . 8) i)
+ (aref a (+ off 7)) (ldb '(8 . 0) i)))
+
+(defun decode-ui64 (a off)
+ "Reads a eight byte integer from the byte array in big endian format."
+ (logior (ash (aref a off) 56)
+ (ash (aref a (+ off 1)) 48)
+ (ash (aref a (+ off 2)) 40)
+ (ash (aref a (+ off 3)) 32)
+ (ash (aref a (+ off 4)) 24)
+ (ash (aref a (+ off 5)) 16)
+ (ash (aref a (+ off 6)) 8)
+ (aref a (+ off 7))))
+
+(defun read-ui32 (stream)
+ (let ((x (make-array 4 :element-type 'unsigned-byte)))
+ (read-sequence x stream)
+ (decode-ui32 x 0)))
+
+;;;; Encoding double float matrices
+
+(define-constant +lisplab-dump-nonce+ 154777230)
+
+(define-constant +lisplab-dump-dge+ 1025)
+
+(defun encode-dge-hdr (rows cols)
+ ;; nonce type skip .... rows cols
+ (let ((x (make-array (* 5 4) :element-type 'unsigned-byte)))
+ (encode-ui32 x +lisplab-dump-nonce+ 0)
+ (encode-ui32 x +lisplab-dump-dge+ 4)
+ (encode-ui32 x 0 8)
+ (encode-ui32 x rows 12)
+ (encode-ui32 x cols 16)
+ x))
+
+(defun encode-dge-bulk (x)
+ (let* ((len (length x))
+ (a (make-array (* 8 len) :element-type 'unsigned-byte)))
+ (dotimes (i len)
+ (encode-ui64 a (ieee-floats:encode-float64 (aref x i)) (* i 8)))
+ a))
+
+(defmethod msave ((s stream) (a matrix-base-dge))
+ ;; Only for binary streams
+ (let ((store (vector-store a)))
+ (write-sequence (encode-dge-hdr (rows a) (cols a)) s)
+ (write-sequence (encode-dge-bulk store) s))
+ (values))
+
+;;; Decoding double float matrices
+
+(defmethod mload ((stream stream))
+ (if (/= (read-ui32 stream) +lisplab-dump-nonce+)
+ (error "Unknown file format for mload.")
+ (if (/= (read-ui32 stream) +lisplab-dump-dge+)
+ (error "Unknown file format for mload.")
+ (progn
+ (let ((len (read-ui32 stream)))
+ (dotimes (i len) (read-byte stream)))
+ (let* ((rows (read-ui32 stream))
+ (cols (read-ui32 stream))
+ (len (* rows cols))
+ (data (make-array (* 8 len) :element-type 'unsigned-byte))
+ (store (allocate-real-store len)))
+ (read-sequence data stream)
+ (dotimes (i len)
+ (setf (aref store i) (ieee-floats:decode-float64
+ (decode-ui64 data (* 8 i)))))
+ (make-instance 'matrix-dge :dim (list rows cols) :store store ))))))
+
+
\ No newline at end of file
More information about the lisplab-cvs
mailing list