diff --git a/CLAUDE.md b/CLAUDE.md index 467dbbd..d5a0b07 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -59,8 +59,9 @@ Tests are split into feature systems, each of which can be run on its own: | `cl-random-forest-test/parallel` | parallelized training accuracy (4, SBCL only) | | `cl-random-forest-test/regression` | univariate regression behaviour (5) | | `cl-random-forest-test/pruning` | global pruning behaviour (5) | +| `cl-random-forest-test/packed` | packed inference representation (11) | -`cl-random-forest-test` is the aggregate that runs all seven. +`cl-random-forest-test` is the aggregate that runs all eight. `cl-random-forest-test/fixture` holds the shared dataset loaders and helpers and has no tests. Load the feature system first (`ql:quickload` or `asdf:load-system`), then: @@ -71,7 +72,13 @@ Load the feature system first (`ql:quickload` or `asdf:load-system`), then: ``` There is no lint step. CI (`.github/workflows/ci.yml`) runs the matrix -{sbcl-bin, ccl-bin} × {ubuntu-latest, macOS-latest}. +{sbcl-bin, ccl-bin} × {ubuntu-latest, macOS-latest}, less ccl-bin on macOS, which the +workflow excludes because macOS CCL binaries are no longer distributed — three jobs, not +four. Roswell has no ccl-bin for **ARM64 Linux** either, so on an aarch64 development +machine the CCL half of the matrix cannot be reproduced locally at all and CI is the only +place it runs. Treat a green local suite as evidence about SBCL only: CCL does not check +type declarations, and its `integer-decode-float` normalises a denormal's significand where +SBCL leaves it alone — both have produced CI-only failures in this repository. Test/example caveats: - `cl-random-forest-test/regression`, `.../pruning`, and the seven synthetic tests in @@ -79,6 +86,7 @@ Test/example caveats: and three `refine-learner-process-*`) need no network: they use `cl-random-forest-test/fixture`'s deterministic synthetic data, or build their own. Everything else downloads datasets. +- `cl-random-forest-test/packed` uses the fixture's synthetic data and needs no network. - Of the regression and pruning suites' ten tests, seven are **property assertions**, not pinned accuracy numbers, and three deliberately pin bugs that are still open: `regression-refine-learner-default-gamma-diverges` (issue #16), `pruning-strands-leaves-without-sample-indices` (issue #14) and @@ -185,6 +193,49 @@ refine learner before training again. the kernel to `nil` restores serial execution with no recompilation. Parallelized: `make-forest`, `make-regression-forest`, `make-refine-dataset`, `train-refine-learner`. +### Packed inference + +`src/packed/` is a separate system, `cl-random-forest/src/packed`, that the core does not +depend on and the facade does not re-export. Load it explicitly: + +```lisp +(ql:quickload :cl-random-forest/src/packed) +``` + +It flattens a trained forest into arrays for inference: 5-8x faster than `predict-forest`, +reentrant where `predict-forest` is not, and serialisable. It is a derived read-only view -- +training, pruning, feature importance and reconstruction all keep using the `node` structs. + +Two things to know. A packed model is a snapshot: prune the forest and the packed copy goes +on predicting with the old structure, silently. And `build-packed-topology` refuses a forest +that has been *pruned*: `pruning!`'s `delete-children!` turns an already-split node -- whose +sample indices were nil'd when it was split -- back into a leaf without restoring them, so +that leaf's class distribution would come out uniform rather than signalling (issue #14). A +default-built, unpruned forest packs fine regardless of what `:remove-sample-indices?` was +given at construction time; `t/packed.lisp` asserts this directly. + +The whole walk runs at `(safety 0)`, so **every index it will follow is checked before it +starts** -- in `check-packable` when building, and again in `packed-load` for a file the +process did not write. Three of those checks are not obvious: + +- `left`/`right` must not merely be in range but must *exceed their own node index*. + `build-packed-topology` numbers a subtree after its root, so this holds by construction; + it is what makes the walk provably terminate, and without it `left[i] = i` is an in-range + cycle that hangs prediction with no way to recover. +- `feature` is bounded by `datum-dim`, the training width, which is a slot on the topology + and a field in the file for exactly this reason. `datum-dim` is itself untrusted in a + file, so it is paired with `check-datamatrix-width`, which bounds it by the matrix + actually passed in. Composed, `feature < datum-dim <= the real column count`. +- Each array's declared length is checked against the bytes left in the file *before* it is + allocated, so an inflated count is an error rather than a multi-gigabyte allocation. + +`check-datamatrix-width` also catches the ordinary version of that mistake -- predicting +with a narrower dataset than the forest was trained on -- which is otherwise a silent read +past the end of a row. It costs nothing measurable: 39-41k predictions/s on a 500-tree +depth-10 forest either way, a spread smaller than the run-to-run noise. + +The design and the measurements behind it are in `docs/packed-forest-layout.md`. + ## Known broken code `src/feature-importance.lisp` and `src/reconstruction.lisp` `:use` only the *exported* symbols of diff --git a/README.org b/README.org index 7f73d08..dcf5ecf 100644 --- a/README.org +++ b/README.org @@ -182,6 +182,32 @@ The following figure shows the accuracy for test dataset and the number of leaf [[./docs/img/clrf-mnist-pruning.png]] +**** Packed inference + +For serving a trained model, =cl-random-forest/src/packed= flattens a forest into arrays: +5-8x faster than =predict-forest=, safe to call from several threads, and serialisable. + +A forest you intend to prune /and/ pack must be built with =:remove-sample-indices? nil=: +pruning turns internal nodes back into leaves without restoring their sample indices (issue +#14), and the packed builder refuses those. =*forest*= above was built with the default +=:remove-sample-indices? t= and was then pruned in the Global Pruning section, so packing it +as shown below will signal a =packed-build-error= unless it is rebuilt with +=:remove-sample-indices? nil= first. + +#+BEGIN_SRC lisp +(ql:quickload :cl-random-forest/src/packed) + +(defparameter *packed* (clrf.packed:build-packed-classifier *forest*)) +(defparameter *acc* (clrf.packed:make-packed-accumulator *packed*)) + +(clrf.packed:packed-predict *packed* *datamatrix* 0 *acc*) + +(clrf.packed:packed-save *packed* #p"model.packed") +(defparameter *reloaded* (clrf.packed:packed-load #p"model.packed")) +#+END_SRC + +One accumulator per thread. The packed model is a snapshot of the forest at build time. + **** Parallelization The following several functions can be parallelized with [[https://lparallel.org/][lparallel]]. diff --git a/cl-random-forest-test.asd b/cl-random-forest-test.asd index 143dbe8..5a81ecd 100644 --- a/cl-random-forest-test.asd +++ b/cl-random-forest-test.asd @@ -65,6 +65,14 @@ (unless (symbol-call :rove :run c) (error "Tests failed.")))) +(defsystem "cl-random-forest-test/packed" + :description "Tests for the packed inference representation" + :depends-on ("rove" "cl-random-forest-test/fixture" "cl-random-forest/src/packed") + :components ((:module "t" :components ((:file "packed")))) + :perform (test-op (o c) (declare (ignore o)) + (unless (symbol-call :rove :run c) + (error "Tests failed.")))) + (defsystem "cl-random-forest-test" :author "Satoshi Imai" :license "MIT Licence" @@ -76,7 +84,8 @@ "cl-random-forest-test/refinement" "cl-random-forest-test/parallel" "cl-random-forest-test/regression" - "cl-random-forest-test/pruning") + "cl-random-forest-test/pruning" + "cl-random-forest-test/packed") ;; NOTE: This system has no components of its own, so rove's SYSTEM-SUITES would ;; return an empty list and (rove:run c) would silently report zero tests. ;; The feature systems must be listed explicitly. @@ -89,5 +98,6 @@ "cl-random-forest-test/refinement" "cl-random-forest-test/parallel" "cl-random-forest-test/regression" - "cl-random-forest-test/pruning")) + "cl-random-forest-test/pruning" + "cl-random-forest-test/packed")) (error "Tests failed.")))) diff --git a/cl-random-forest.asd b/cl-random-forest.asd index 3461ff3..15c4f9b 100644 --- a/cl-random-forest.asd +++ b/cl-random-forest.asd @@ -2,6 +2,18 @@ (:use :cl :asdf)) (in-package :cl-random-forest-asd) +;; SVMFORMAT is a nickname of the CL-LIBSVM-FORMAT package, and package-inferred-system +;; derives a system name by downcasing the package name -- which would look for a +;; nonexistent "svmformat" system. Registering the mapping is what lets src/utils.lisp +;; declare its dependency on the package it actually uses. +(asdf:register-system-packages "cl-libsvm-format" '(#:svmformat)) + +;; CLOL and CLOL.VECTOR are nicknames of CL-ONLINE-LEARNING and CL-ONLINE-LEARNING.VECTOR +;; respectively, both provided by the single CL-ONLINE-LEARNING system. Same reasoning as +;; the CL-LIBSVM-FORMAT mapping above: without this, src/random-forest.lisp's :import-from +;; of these nicknames would send ASDF looking for nonexistent "clol"/"clol.vector" systems. +(asdf:register-system-packages "cl-online-learning" '(#:clol #:clol.vector)) + (defsystem cl-random-forest :version "0.2" :author "Satoshi Imai" diff --git a/src/packed.lisp b/src/packed.lisp new file mode 100644 index 0000000..c0db616 --- /dev/null +++ b/src/packed.lisp @@ -0,0 +1,8 @@ +(uiop:define-package :cl-random-forest/src/packed + (:use :cl) + (:nicknames :clrf.packed) + (:use-reexport :cl-random-forest/src/packed/topology) + (:use-reexport :cl-random-forest/src/packed/classifier) + (:use-reexport :cl-random-forest/src/packed/io)) + +(in-package :cl-random-forest/src/packed) diff --git a/src/packed/classifier.lisp b/src/packed/classifier.lisp new file mode 100644 index 0000000..6068de3 --- /dev/null +++ b/src/packed/classifier.lisp @@ -0,0 +1,266 @@ +(defpackage :cl-random-forest/src/packed/classifier + (:use #:cl + #:cl-random-forest/src/packed/topology) + (:import-from #:cl-random-forest/src/random-forest + #:forest-dtree-list + #:forest-n-class + #:dtree-root + #:node-test-attribute + #:node-left-node + #:node-right-node + #:node-class-distribution + #:class-distribution-forest + #:predict-forest + #:argmax) + (:export #:packed-classifier + #:packed-classifier-p + #:packed-classifier-topology + #:packed-classifier-n-class + #:packed-classifier-kind + #:packed-classifier-table + #:packed-classifier-offsets + #:packed-classifier-class + #:packed-classifier-probability + #:%make-packed-classifier + #:build-packed-classifier + #:make-packed-accumulator + #:make-packed-accumulators + #:packed-predict + #:packed-predict-batch + #:packed-verify)) + +(in-package :cl-random-forest/src/packed/classifier) + +(defstruct (packed-classifier (:constructor %make-packed-classifier)) + "A topology plus a leaf payload. + +KIND is :DENSE, in which case TABLE is an n-leaf x n-class array and the CSR slots are +empty; or :CSR, in which case leaf L occupies [OFFSETS[L], OFFSETS[L+1]) of CLASS and +PROBABILITY and TABLE is empty. Omitting the zeros is exact, not approximate: the values +are non-negative, so leaving out a + 0.0 changes no sum." + (topology (%make-packed-topology) :type packed-topology) + (n-class 0 :type fixnum) + (kind :dense :type (member :dense :csr)) + (table (make-array '(0 0) :element-type 'single-float) + :type (simple-array single-float (* *))) + (offsets (make-array 1 :element-type '(unsigned-byte 32)) + :type (simple-array (unsigned-byte 32) (*))) + (class (make-array 0 :element-type '(unsigned-byte 16)) + :type (simple-array (unsigned-byte 16) (*))) + (probability (make-array 0 :element-type 'single-float) + :type (simple-array single-float (*)))) + +(defun leaf-distributions (forest topology) + "A simple-vector of every leaf's distribution, indexed by the topology's leaf number. + +NODE-CLASS-DISTRIBUTION hands back the dtree's shared scratch buffer, so each one is copied +before the next leaf overwrites it." + (let ((out (make-array (packed-topology-n-leaf topology))) + (offsets (packed-topology-tree-leaf-offsets topology))) + (loop for dtree in (forest-dtree-list forest) + for tree from 0 + do (let ((local 0) + (base (aref offsets tree))) + (labels ((walk (node) + (if (node-test-attribute node) + (progn (walk (node-left-node node)) + (walk (node-right-node node))) + (progn (setf (svref out (+ base local)) + (copy-seq (node-class-distribution node))) + (incf local))))) + (walk (dtree-root dtree))))) + out)) + +(defun choose-payload (n-class) + "CSR once a dense row exceeds a cache line, dense otherwise. + +The tempting rule is the compression ratio, and it picks wrong: a 10-class model compresses +2.4x and runs 0.82x as fast, because a 40-byte dense row already fits in a cache line and +CSR's offset loads and scattered writes cost more than the bytes saved." + (if (> (* 4 n-class) 64) :csr :dense)) + +(defun build-packed-classifier (forest &key (payload :auto)) + "Flatten FOREST into a packed classifier. PAYLOAD is :AUTO, :DENSE or :CSR." + (check-type payload (member :auto :dense :csr)) + (let* ((topology (build-packed-topology forest)) + (n-class (forest-n-class forest)) + (kind (if (eq payload :auto) (choose-payload n-class) payload)) + (n-leaf (packed-topology-n-leaf topology)) + (dists (leaf-distributions forest topology))) + (ecase kind + (:dense + (let ((table (make-array (list (max n-leaf 1) n-class) + :element-type 'single-float :initial-element 0.0))) + (dotimes (row n-leaf) + (let ((d (svref dists row))) + (dotimes (k n-class) (setf (aref table row k) (aref d k))))) + (%make-packed-classifier :topology topology :n-class n-class + :kind :dense :table table))) + (:csr + (unless (< n-class 65536) + (error 'packed-build-error + :detail (format nil "~D classes exceeds a (unsigned-byte 16) class index" + n-class))) + (let ((nnz 0)) + (dotimes (row n-leaf) + (let ((d (svref dists row))) + (dotimes (k n-class) (unless (zerop (aref d k)) (incf nnz))))) + (let ((offsets (make-array (1+ n-leaf) :element-type '(unsigned-byte 32))) + (class (make-array (max nnz 1) :element-type '(unsigned-byte 16))) + (probability (make-array (max nnz 1) :element-type 'single-float)) + (cursor 0)) + (dotimes (row n-leaf) + (setf (aref offsets row) cursor) + (let ((d (svref dists row))) + (dotimes (k n-class) + (let ((p (aref d k))) + (unless (zerop p) + (setf (aref class cursor) k + (aref probability cursor) p) + (incf cursor)))))) + (setf (aref offsets n-leaf) cursor) + (%make-packed-classifier :topology topology :n-class n-class :kind :csr + :offsets offsets :class class + :probability probability))))))) + +(defun make-packed-accumulator (classifier) + "A fresh accumulator for PACKED-PREDICT. One per thread, never shared." + (make-array (packed-classifier-n-class classifier) + :element-type 'single-float :initial-element 0.0)) + +(defun packed-predict (classifier datamatrix datum-index acc) + "PREDICT-FOREST's answer for one datum. Writes only ACC. + +ACC is normalised in place: on return it holds the class distribution (summed votes divided +by N-TREE), not the raw sums. PACKED-VERIFY relies on this to read the distribution straight +out of ACC. Contrast PACKED-PREDICT-BATCH, which leaves its accumulators unnormalised." + (declare (optimize (speed 3) (safety 0)) + (type packed-classifier classifier) + (type (simple-array single-float (* *)) datamatrix) + (type (simple-array single-float (*)) acc) + (type fixnum datum-index)) + (let* ((topology (packed-classifier-topology classifier)) + (n-class (packed-classifier-n-class classifier)) + (roots (packed-topology-roots topology)) + (n-tree (packed-topology-n-tree topology))) + (declare (type (simple-array (signed-byte 32) (*)) roots) + (type fixnum n-class n-tree)) + (check-datamatrix-width topology datamatrix) + (dotimes (k n-class) (setf (aref acc k) 0.0)) + (ecase (packed-classifier-kind classifier) + (:dense + (let ((table (packed-classifier-table classifier))) + (declare (type (simple-array single-float (* *)) table)) + (dotimes (tree n-tree) + (let ((row (packed-leaf topology datamatrix datum-index (aref roots tree)))) + (declare (type fixnum row)) + (dotimes (k n-class) (incf (aref acc k) (aref table row k))))))) + (:csr + (let ((offsets (packed-classifier-offsets classifier)) + (class (packed-classifier-class classifier)) + (probability (packed-classifier-probability classifier))) + (declare (type (simple-array (unsigned-byte 32) (*)) offsets) + (type (simple-array (unsigned-byte 16) (*)) class) + (type (simple-array single-float (*)) probability)) + (dotimes (tree n-tree) + (let ((leaf (packed-leaf topology datamatrix datum-index (aref roots tree)))) + (declare (type fixnum leaf)) + (loop for i of-type fixnum + from (aref offsets leaf) below (aref offsets (1+ leaf)) + do (incf (aref acc (aref class i)) (aref probability i)))))))) + (dotimes (k n-class) (setf (aref acc k) (/ (aref acc k) n-tree))) + (argmax acc))) + +(defun make-packed-accumulators (classifier tile) + "Accumulators for PACKED-PREDICT-BATCH: one row per datum in the tile." + (make-array (list tile (packed-classifier-n-class classifier)) + :element-type 'single-float :initial-element 0.0)) + +(defun packed-predict-batch (classifier datamatrix start end accs out) + "Predict rows [START,END) a tree at a time, writing classes into OUT indexed from zero. + +Walking one tree over the whole tile before moving to the next touches each tree's arrays +once per tile instead of once per datum. The cost is holding (END - START) x n-class +accumulators, which is why a tile of one loses and a tile of a few hundred wins. + +ACCS must have at least (END - START) rows -- as MAKE-PACKED-ACCUMULATORS gives it when +called with a TILE no smaller than (END - START) -- and OUT must have length at least +(END - START); this runs at (safety 0), so either being too small is a silent out-of-bounds +write, not a signalled error. + +Unlike PACKED-PREDICT, ACCS is left holding raw per-tree sums, not divided by N-TREE: the +division by N-TREE happens only in a local while picking each row's argmax for OUT, and +that quotient is never written back. A caller that wants ACCS's own rows to be class +distributions must divide each one by N-TREE itself." + (declare (optimize (speed 3) (safety 0)) + (type packed-classifier classifier) + (type (simple-array single-float (* *)) datamatrix accs) + (type (simple-array fixnum (*)) out) + (type fixnum start end)) + (let* ((topology (packed-classifier-topology classifier)) + (n-class (packed-classifier-n-class classifier)) + (roots (packed-topology-roots topology)) + (n-tree (packed-topology-n-tree topology)) + (tile (- end start))) + (declare (type (simple-array (signed-byte 32) (*)) roots) + (type fixnum n-class n-tree tile)) + (check-datamatrix-width topology datamatrix) + (dotimes (r tile) + (dotimes (k n-class) (setf (aref accs r k) 0.0))) + (ecase (packed-classifier-kind classifier) + (:dense + (let ((table (packed-classifier-table classifier))) + (declare (type (simple-array single-float (* *)) table)) + (dotimes (tree n-tree) + (let ((root (aref roots tree))) + (dotimes (r tile) + (let ((row (packed-leaf topology datamatrix (+ start r) root))) + (declare (type fixnum row)) + (dotimes (k n-class) + (incf (aref accs r k) (aref table row k))))))))) + (:csr + (let ((offsets (packed-classifier-offsets classifier)) + (class (packed-classifier-class classifier)) + (probability (packed-classifier-probability classifier))) + (declare (type (simple-array (unsigned-byte 32) (*)) offsets) + (type (simple-array (unsigned-byte 16) (*)) class) + (type (simple-array single-float (*)) probability)) + (dotimes (tree n-tree) + (let ((root (aref roots tree))) + (dotimes (r tile) + (let ((leaf (packed-leaf topology datamatrix (+ start r) root))) + (declare (type fixnum leaf)) + (loop for i of-type fixnum + from (aref offsets leaf) below (aref offsets (1+ leaf)) + do (incf (aref accs r (aref class i)) + (aref probability i)))))))))) + (dotimes (r tile out) + (let ((best most-negative-single-float) + (best-k 0)) + (declare (type single-float best) (type fixnum best-k)) + (dotimes (k n-class) + (let ((v (/ (aref accs r k) n-tree))) + (when (> v best) (setf best v best-k k)))) + (setf (aref out r) best-k))))) + +(defun packed-verify (classifier forest datamatrix) + "Compare CLASSIFIER with FOREST over every row of DATAMATRIX. + +Returns (values class-disagreements distribution-disagreements worst-absolute-difference). +All three should be zero: the distribution is checked element by element, not merely after +ARGMAX, and exact equality is achievable because the same floats are summed in the same +order." + (let* ((n-class (packed-classifier-n-class classifier)) + (acc (make-packed-accumulator classifier)) + (class-bad 0) (dist-bad 0) (worst 0.0)) + (dotimes (i (array-dimension datamatrix 0) (values class-bad dist-bad worst)) + (let ((reference-class (predict-forest forest datamatrix i)) + (reference-dist (copy-seq (class-distribution-forest forest datamatrix i)))) + (unless (= reference-class (packed-predict classifier datamatrix i acc)) + (incf class-bad)) + (let ((row-bad nil)) + (dotimes (k n-class) + (let ((d (abs (- (aref reference-dist k) (aref acc k))))) + (when (> d 0.0) (setf row-bad t)) + (setf worst (max worst d)))) + (when row-bad (incf dist-bad))))))) diff --git a/src/packed/io.lisp b/src/packed/io.lisp new file mode 100644 index 0000000..a454717 --- /dev/null +++ b/src/packed/io.lisp @@ -0,0 +1,406 @@ +(defpackage :cl-random-forest/src/packed/io + (:use #:cl + #:cl-random-forest/src/packed/topology + #:cl-random-forest/src/packed/classifier) + (:export #:packed-save + #:packed-load + #:packed-load-error + #:packed-load-error-detail + #:single-float-to-bits + #:bits-to-single-float + #:%portable-single-float-to-bits + #:%portable-bits-to-single-float + #:%float-parts-to-bits)) + +(in-package :cl-random-forest/src/packed/io) + +(define-condition packed-load-error (error) + ((detail :initarg :detail :reader packed-load-error-detail)) + (:report (lambda (c s) + (format s "cannot load this packed model: ~A" (packed-load-error-detail c))))) + +;; A DEFCONSTANT whose value is an array signals on reload, the new value not being EQL to +;; the old, so the magic is a DEFPARAMETER. The other two are numbers and are fine. +(defparameter *magic* + #.(map '(simple-array (unsigned-byte 8) (*)) #'char-code "CLRFPACK")) +(defconstant +version+ 2) +(defconstant +byte-order-probe+ #x01020304) + +;;;; Floats as bits +;;;; +;;;; The file holds IEEE-754 bit patterns. SBCL and CCL can produce them directly; the +;;;; portable version is the reference the tests hold the fast paths to, and the fallback +;;;; anywhere else. + +(defun %float-parts-to-bits (significand exponent sign) + "IEEE-754 single-precision bits of SIGNIFICAND * 2^EXPONENT, signed by SIGN. + +Separate from %PORTABLE-SINGLE-FLOAT-TO-BITS so that both of the conventions +INTEGER-DECODE-FLOAT is allowed to use can be fed in from any host and checked, rather than +only on a host that happens to use the other one." + (let ((sign-bit (if (minusp sign) #x80000000 0))) + ;; INTEGER-DECODE-FLOAT owes its caller only a pair whose product is the magnitude; + ;; whether a denormal's significand comes back normalised is up to the implementation. + ;; SBCL leaves it alone -- LEAST-POSITIVE-SINGLE-FLOAT is significand 1, exponent -149 -- + ;; and CCL normalises it to significand 2^23, exponent -172. Deciding "is this a + ;; denormal?" by the significand's magnitude therefore reads the host's convention + ;; instead of the format's, which is exactly how this went wrong before: correct on SBCL, + ;; and on CCL it fell into the normal-number branch and wrote a wrapped exponent field. + ;; Normalising first and letting the biased exponent decide is true of either convention. + ;; The pair is pinned down only up to doubling the significand and dropping the exponent, + ;; so reduce it from both directions rather than assuming which end it arrived at. + (loop while (and (>= significand (ash 1 24)) (evenp significand)) + do (setf significand (ash significand -1)) + (incf exponent)) + (loop while (< significand (ash 1 23)) + do (setf significand (ash significand 1)) + (decf exponent)) + (let ((biased (+ exponent 23 127))) + (if (plusp biased) + (logior sign-bit (ash biased 23) (logand significand #x7fffff)) + ;; Exponent field zero, and the significand shifted down to the fixed 2^-149 scale + ;; a denormal's fraction is measured in. + (logior sign-bit (ash significand (1- biased))))))) + +(defun %portable-single-float-to-bits (x) + "IEEE-754 single-precision bits of X, in ANSI Common Lisp only." + (declare (type single-float x)) + (if (zerop x) + (if (minusp (float-sign x)) #x80000000 0) + (multiple-value-bind (significand exponent sign) (integer-decode-float x) + (%float-parts-to-bits significand exponent sign)))) + +(defun %portable-bits-to-single-float (bits) + "The single-float whose IEEE-754 bits are BITS, in ANSI Common Lisp only." + (declare (type (unsigned-byte 32) bits)) + (let* ((sign (if (logbitp 31 bits) -1.0 1.0)) + (exponent (ldb (byte 8 23) bits)) + (fraction (ldb (byte 23 0) bits))) + (cond ((and (zerop exponent) (zerop fraction)) (* sign 0.0)) + ((zerop exponent) (* sign (scale-float (float fraction 1.0) -149))) + (t (* sign (scale-float (float (logior fraction (ash 1 23)) 1.0) + (- exponent 127 23))))))) + +(defun single-float-to-bits (x) + "IEEE-754 single-precision bits of X." + (declare (type single-float x)) + #+sbcl (ldb (byte 32 0) (sb-kernel:single-float-bits x)) + #+ccl (ldb (byte 32 0) (ccl::single-float-bits x)) + #-(or sbcl ccl) (%portable-single-float-to-bits x)) + +(defun bits-to-single-float (bits) + "The single-float whose IEEE-754 bits are BITS." + (declare (type (unsigned-byte 32) bits)) + #+sbcl (sb-kernel:make-single-float + (if (logbitp 31 bits) (- bits (ash 1 32)) bits)) + #+ccl (ccl::host-single-float-from-unsigned-byte-32 bits) + #-(or sbcl ccl) (%portable-bits-to-single-float bits)) + +;;;; Byte-level primitives + +(defun write-u32 (value stream) + (declare (type (unsigned-byte 32) value)) + (dotimes (i 4) (write-byte (ldb (byte 8 (* 8 i)) value) stream))) + +(defun read-u32 (stream) + (let ((value 0)) + (dotimes (i 4 value) + (let ((byte (read-byte stream nil nil))) + (unless byte (error 'packed-load-error :detail "the file ends mid-header")) + (setf value (logior value (ash byte (* 8 i)))))))) + +(defun write-array (array stream bits encoder) + "Write ARRAY as little-endian words of BITS bits, ENCODER mapping element to integer." + (let* ((bytes-per (floor bits 8)) + (n (length array)) + (buffer (make-array (* n bytes-per) :element-type '(unsigned-byte 8)))) + (dotimes (i n) + (let ((word (funcall encoder (aref array i)))) + (dotimes (b bytes-per) + (setf (aref buffer (+ (* i bytes-per) b)) (ldb (byte 8 (* 8 b)) word))))) + (write-sequence buffer stream))) + +(defun check-array-fits (stream name n bits) + "Signal PACKED-LOAD-ERROR unless N words of BITS bits still remain in STREAM. + +READ-ARRAY has to size two buffers from N before it can discover that the file is shorter +than N claims, and N is a header field a corrupt file controls outright: a 36-byte file +declaring 4294967295 internal nodes would ask for 16 GB before reaching the short-read check. +Comparing against the bytes actually left costs one FILE-LENGTH and forecloses that." + (let* ((wanted (* n (floor bits 8))) + (remaining (- (file-length stream) (file-position stream)))) + (when (> wanted remaining) + (error 'packed-load-error + :detail (format nil "~A claims ~D entries, needing ~D bytes, but only ~D ~ +remain in the file" name n wanted remaining))))) + +(defun read-array (stream n bits element-type decoder) + "Read N little-endian words of BITS bits, DECODER mapping integer to element." + (let* ((bytes-per (floor bits 8)) + (buffer (make-array (* n bytes-per) :element-type '(unsigned-byte 8))) + (got (read-sequence buffer stream)) + (out (make-array n :element-type element-type))) + (unless (= got (* n bytes-per)) + (error 'packed-load-error + :detail (format nil "wanted ~D bytes of array data, got ~D" + (* n bytes-per) got))) + (dotimes (i n out) + (let ((word 0)) + (dotimes (b bytes-per) + (setf word (logior word (ash (aref buffer (+ (* i bytes-per) b)) (* 8 b))))) + (setf (aref out i) (funcall decoder word)))))) + +(defun u32-of-s32 (x) (ldb (byte 32 0) x)) +(defun s32-of-u32 (x) (if (logbitp 31 x) (- x (ash 1 32)) x)) + +;;;; Loaded-array validation +;;;; +;;;; BUILD-PACKED-TOPOLOGY validates exhaustively before anything reaches its (safety 0) +;;;; traversal. PACKED-LOAD reconstructs the same arrays from bytes it does not control, so +;;;; it owes them the same discipline: every index that PACKED-LEAF or the CSR loop in +;;;; PACKED-PREDICT will follow at (safety 0) is checked here first, where a bad one is an +;;;; ordinary error instead of an out-of-bounds read or write. Thresholds and the leaf +;;;; payload are not checked -- any bit pattern is a legal float, and NaN is not this +;;;; layer's problem. + +(defun check-header-count (name value &key positive) + "Signal PACKED-LOAD-ERROR unless VALUE is a valid header count: non-negative always, and +positive too when POSITIVE is true. READ-U32 cannot actually hand back a negative value, but +the check is stated explicitly rather than left as an accident of the encoding." + (unless (<= 0 value) + (error 'packed-load-error :detail (format nil "header ~A is negative: ~D" name value))) + (when (and positive (not (plusp value))) + (error 'packed-load-error + :detail (format nil "header ~A must be positive, got ~D" name value)))) + +(defun check-child-reference (array-name index value n-internal n-leaf &key parent) + "Signal PACKED-LOAD-ERROR unless VALUE is a valid child reference at ARRAY-NAME[INDEX]: +either an internal-node index in [0, N-INTERNAL), or the BUILD-PACKED-TOPOLOGY encoding of a +leaf, `~L` for some leaf L in [0, N-LEAF). + +When PARENT is the index of the node VALUE is a child of, an internal VALUE must also exceed +it. Bounds alone do not make the walk terminate: PACKED-LEAF loops while the node index is +non-negative, so `left[i] = i` -- in range, and a cycle -- hangs it at (safety 0). Requiring +each step to move forward through a bounded array makes termination provable instead, which +is exactly the numbering BUILD-PACKED-TOPOLOGY produces and now checks." + (if (minusp value) + (let ((leaf (lognot value))) + (unless (< -1 leaf n-leaf) + (error 'packed-load-error + :detail (format nil "~A[~D] = ~D references leaf ~D, outside [0,~D)" + array-name index value leaf n-leaf)))) + (progn + (unless (< value n-internal) + (error 'packed-load-error + :detail (format nil "~A[~D] = ~D, outside the internal-node range [0,~D)" + array-name index value n-internal))) + (when (and parent (<= value parent)) + (error 'packed-load-error + :detail (format nil "~A[~D] = ~D does not exceed its parent index ~D, so ~ +the walk could revisit a node and never terminate" array-name index value parent)))))) + +(defun validate-feature-array (n-internal datum-dim feature) + "Signal PACKED-LOAD-ERROR unless every FEATURE entry is a column of a DATUM-DIM datamatrix. + +PACKED-LEAF uses FEATURE[node] as the second subscript of the datamatrix with bounds checking +off, so an unchecked entry from a corrupt file is a read at an arbitrary offset from the +array. CHECK-PACKABLE applies the same bound when building; the file carries DATUM-DIM so +that the bound survives the round trip. + +DATUM-DIM is itself a number from the file, so this check alone would be vacuous against a +corruption that inflated both. It is half of a pair: this bounds FEATURE by the claimed +width, and CHECK-DATAMATRIX-WIDTH bounds the claimed width by the matrix actually handed to +the walk. Composed, feature < datum-dim <= the real column count, and an inflated DATUM-DIM +fails closed at the first prediction rather than reading out of bounds." + (dotimes (i n-internal) + (unless (< (aref feature i) datum-dim) + (error 'packed-load-error + :detail (format nil "feature[~D] = ~D, outside the datamatrix's [0,~D) columns" + i (aref feature i) datum-dim))))) + +(defun validate-topology-arrays (n-internal n-leaf n-tree left right roots tree-leaf-offsets) + "Signal PACKED-LOAD-ERROR if LEFT, RIGHT, ROOTS or TREE-LEAF-OFFSETS are not internally +consistent with the header counts N-INTERNAL, N-LEAF and N-TREE." + (dotimes (i n-internal) + (check-child-reference "left" i (aref left i) n-internal n-leaf :parent i) + (check-child-reference "right" i (aref right i) n-internal n-leaf :parent i)) + ;; A root has no parent to exceed; it only has to be in range. Termination still follows, + ;; since every step after it moves strictly forward. + (dotimes (tree n-tree) + (check-child-reference "roots" tree (aref roots tree) n-internal n-leaf)) + (unless (zerop (aref tree-leaf-offsets 0)) + (error 'packed-load-error + :detail (format nil "tree-leaf-offsets[0] = ~D, must start at 0" + (aref tree-leaf-offsets 0)))) + (let ((previous 0)) + (dotimes (tree n-tree) + (let ((offset (aref tree-leaf-offsets tree))) + (when (< offset previous) + (error 'packed-load-error + :detail (format nil "tree-leaf-offsets[~D] = ~D is less than the previous ~ +entry ~D" tree offset previous))) + (unless (or (< offset n-leaf) (and (zerop n-leaf) (zerop offset))) + (error 'packed-load-error + :detail (format nil "tree-leaf-offsets[~D] = ~D, outside [0,~D)" + tree offset n-leaf))) + (setf previous offset))))) + +(defun validate-csr-arrays (n-leaf n-class offsets class probability) + "Signal PACKED-LOAD-ERROR unless OFFSETS, CLASS and PROBABILITY are a valid CSR partition +of N-LEAF rows over N-CLASS columns." + (unless (zerop (aref offsets 0)) + (error 'packed-load-error + :detail (format nil "csr offsets[0] = ~D, must start at 0" (aref offsets 0)))) + (let ((previous (aref offsets 0))) + (loop for i from 1 to n-leaf + do (let ((offset (aref offsets i))) + (when (< offset previous) + (error 'packed-load-error + :detail (format nil "csr offsets[~D] = ~D is less than the previous ~ +entry ~D" i offset previous))) + (setf previous offset)))) + (let ((nnz (aref offsets n-leaf))) + (unless (= nnz (length class)) + (error 'packed-load-error + :detail (format nil "csr offsets[~D] = ~D but class has ~D entries" + n-leaf nnz (length class)))) + (unless (= nnz (length probability)) + (error 'packed-load-error + :detail (format nil "csr offsets[~D] = ~D but probability has ~D entries" + n-leaf nnz (length probability))))) + (dotimes (i (length class)) + (unless (< (aref class i) n-class) + (error 'packed-load-error + :detail (format nil "csr class[~D] = ~D, outside [0,~D)" + i (aref class i) n-class))))) + +;;;; Save and load + +(defun packed-save (classifier pathname) + "Write CLASSIFIER to PATHNAME. Returns PATHNAME. + +The format is a header followed by the arrays as little-endian words: WRITE-U32 always +serialises LSB-first and SINGLE-FLOAT-TO-BITS yields an endianness-independent IEEE-754 bit +pattern, so the file is byte-order independent by construction and round-trips identically +on any machine. +BYTE-ORDER-PROBE+ is not there to catch a real mismatch, then -- it is a +second, distinctive magic number, cheap insurance that the header was not merely truncated +or shifted past the point the first magic check looks at." + (let* ((topology (packed-classifier-topology classifier)) + (n-leaf (packed-topology-n-leaf topology)) + (n-class (packed-classifier-n-class classifier)) + (csr (eq (packed-classifier-kind classifier) :csr))) + (with-open-file (s pathname :direction :output :if-exists :supersede + :element-type '(unsigned-byte 8)) + (write-sequence *magic* s) + (write-u32 +version+ s) + (write-u32 +byte-order-probe+ s) + (write-u32 (if csr 1 0) s) + (write-u32 (packed-topology-n-tree topology) s) + (write-u32 (packed-topology-n-internal topology) s) + (write-u32 n-leaf s) + (write-u32 n-class s) + (write-u32 (packed-topology-datum-dim topology) s) + (write-array (packed-topology-feature topology) s 32 #'identity) + (write-array (packed-topology-threshold topology) s 32 #'single-float-to-bits) + (write-array (packed-topology-left topology) s 32 #'u32-of-s32) + (write-array (packed-topology-right topology) s 32 #'u32-of-s32) + (write-array (packed-topology-roots topology) s 32 #'u32-of-s32) + (write-array (packed-topology-tree-leaf-offsets topology) s 32 #'identity) + (if csr + (progn + (write-array (packed-classifier-offsets classifier) s 32 #'identity) + (write-array (packed-classifier-class classifier) s 16 #'identity) + (write-array (packed-classifier-probability classifier) s 32 + #'single-float-to-bits)) + (let* ((table (packed-classifier-table classifier)) + (flat (make-array (* (max n-leaf 1) n-class) + :element-type 'single-float))) + (dotimes (row (max n-leaf 1)) + (dotimes (k n-class) + (setf (aref flat (+ (* row n-class) k)) (aref table row k)))) + (write-array flat s 32 #'single-float-to-bits)))) + pathname)) + +(defun packed-load (pathname) + "Read a packed classifier written by PACKED-SAVE. + +Every array is checked against the header counts before this returns: FEATURE against the +saved datamatrix width, LEFT, RIGHT and ROOTS for range and for forward motion, plus +TREE-LEAF-OFFSETS and, for a CSR file, OFFSETS and CLASS. A corrupt file the topology +traversal or the CSR prediction loop would otherwise read or write out of bounds at +(safety 0) -- or walk in a cycle forever -- instead signals PACKED-LOAD-ERROR here. Each +array's declared length is also checked against the bytes actually remaining before it is +allocated, so an inflated count cannot turn into a huge allocation." + (with-open-file (s pathname :direction :input :element-type '(unsigned-byte 8)) + (let ((magic (make-array (length *magic*) :element-type '(unsigned-byte 8)))) + (unless (and (= (read-sequence magic s) (length *magic*)) + (equalp magic *magic*)) + (error 'packed-load-error :detail "not a packed forest file"))) + (let ((version (read-u32 s))) + (unless (= version +version+) + (error 'packed-load-error + :detail (format nil "version ~D, this build reads ~D" version +version+)))) + (let ((order (read-u32 s))) + (unless (= order +byte-order-probe+) + (error 'packed-load-error + :detail (format nil "written on a machine of different byte order (~8,'0X)" + order)))) + (let* ((csr (= 1 (read-u32 s))) + (n-tree (read-u32 s)) + (n-internal (read-u32 s)) + (n-leaf (read-u32 s)) + (n-class (read-u32 s)) + (datum-dim (read-u32 s))) + (check-header-count "n-tree" n-tree :positive t) + (check-header-count "n-internal" n-internal) + (check-header-count "n-leaf" n-leaf) + (check-header-count "n-class" n-class :positive t) + (check-header-count "datum-dim" datum-dim :positive t) + (check-array-fits s "feature" n-internal 32) + (let* ((feature (read-array s n-internal 32 '(unsigned-byte 32) #'identity)) + (threshold (progn (check-array-fits s "threshold" n-internal 32) + (read-array s n-internal 32 'single-float + #'bits-to-single-float))) + (left (progn (check-array-fits s "left" n-internal 32) + (read-array s n-internal 32 '(signed-byte 32) #'s32-of-u32))) + (right (progn (check-array-fits s "right" n-internal 32) + (read-array s n-internal 32 '(signed-byte 32) #'s32-of-u32))) + (roots (progn (check-array-fits s "roots" n-tree 32) + (read-array s n-tree 32 '(signed-byte 32) #'s32-of-u32))) + (tree-leaf-offsets (progn (check-array-fits s "tree-leaf-offsets" n-tree 32) + (read-array s n-tree 32 '(unsigned-byte 32) + #'identity)))) + (validate-topology-arrays n-internal n-leaf n-tree left right roots + tree-leaf-offsets) + (validate-feature-array n-internal datum-dim feature) + (let ((topology + (%make-packed-topology + :n-tree n-tree :n-internal n-internal :n-leaf n-leaf + :datum-dim datum-dim + :feature feature :threshold threshold :left left :right right + :roots roots :tree-leaf-offsets tree-leaf-offsets))) + (if csr + (let* ((offsets (progn (check-array-fits s "csr offsets" (1+ n-leaf) 32) + (read-array s (1+ n-leaf) 32 '(unsigned-byte 32) + #'identity))) + (nnz (aref offsets n-leaf)) + (class (progn (check-array-fits s "csr class" nnz 16) + (read-array s nnz 16 '(unsigned-byte 16) #'identity))) + (probability (progn (check-array-fits s "csr probability" nnz 32) + (read-array s nnz 32 'single-float + #'bits-to-single-float)))) + (validate-csr-arrays n-leaf n-class offsets class probability) + (%make-packed-classifier + :topology topology :n-class n-class :kind :csr + :offsets offsets :class class :probability probability)) + (let ((flat (progn + (check-array-fits s "dense table" (* (max n-leaf 1) n-class) 32) + (read-array s (* (max n-leaf 1) n-class) 32 'single-float + #'bits-to-single-float))) + (table (make-array (list (max n-leaf 1) n-class) + :element-type 'single-float))) + (dotimes (row (max n-leaf 1)) + (dotimes (k n-class) + (setf (aref table row k) (aref flat (+ (* row n-class) k))))) + (%make-packed-classifier + :topology topology :n-class n-class :kind :dense :table table)))))))) diff --git a/src/packed/topology.lisp b/src/packed/topology.lisp new file mode 100644 index 0000000..acdc106 --- /dev/null +++ b/src/packed/topology.lisp @@ -0,0 +1,265 @@ +(defpackage :cl-random-forest/src/packed/topology + (:use #:cl) + (:import-from #:cl-random-forest/src/random-forest + #:forest-dtree-list + #:forest-n-tree + #:forest-datum-dim + #:dtree-root + #:node-test-attribute + #:node-test-threshold + #:node-left-node + #:node-right-node + #:node-sample-indices) + (:export #:packed-topology + #:packed-topology-p + #:packed-topology-n-tree + #:packed-topology-n-internal + #:packed-topology-n-leaf + #:packed-topology-datum-dim + #:packed-topology-feature + #:packed-topology-threshold + #:packed-topology-left + #:packed-topology-right + #:packed-topology-roots + #:packed-topology-tree-leaf-offsets + #:%make-packed-topology + #:build-packed-topology + #:packed-leaf + #:packed-leaf-indices + #:check-datamatrix-width + #:packed-build-error + #:packed-build-error-detail)) + +(in-package :cl-random-forest/src/packed/topology) + +(define-condition packed-build-error (error) + ((detail :initarg :detail :reader packed-build-error-detail)) + (:report (lambda (c s) + (format s "cannot pack this forest: ~A" (packed-build-error-detail c)))) + (:documentation + "Signalled for anything the packed traversal could not survive. + +The traversal runs at (safety 0) and trusts its declarations, so every assumption has to +hold before it starts. One of them is not obvious: a leaf that has lost its sample indices +would get a uniform class distribution from CLASS-DISTRIBUTION rather than an error, and +that would be frozen into the model silently.")) + +(defstruct (packed-topology (:constructor %make-packed-topology)) + "Every internal node of every tree of a forest, flattened into parallel arrays. + +A leaf is encoded as a negative child index, `~leaf`, so a walk's continuation test is its +leaf test and a leaf costs no extra read. Leaf numbers are +TREE-LEAF-OFFSETS[tree] + the tree's own leaf index, which is by construction the index +Global Refinement uses. + +LEFT and RIGHT hold, for every internal node, a child index strictly greater than the node's +own index: EMIT takes its slot before recursing, so a subtree is always numbered after its +root. Nothing in the walk depends on that, but it makes termination structural rather than +circumstantial -- each step moves strictly forward through a bounded array -- and PACKED-LOAD +rechecks it on a file it did not write, where a cycle would otherwise hang the walk. + +DATUM-DIM is the width of the datamatrix the forest was trained on. FEATURE indexes a +datamatrix column at (safety 0), so it is bounded by DATUM-DIM at build time, and keeping the +bound is what lets PACKED-LOAD and the prediction entry points check it again later." + (n-tree 0 :type fixnum) + (n-internal 0 :type fixnum) + (n-leaf 0 :type fixnum) + (datum-dim 0 :type fixnum) + (feature (make-array 0 :element-type '(unsigned-byte 32)) + :type (simple-array (unsigned-byte 32) (*))) + (threshold (make-array 0 :element-type 'single-float) + :type (simple-array single-float (*))) + (left (make-array 0 :element-type '(signed-byte 32)) + :type (simple-array (signed-byte 32) (*))) + (right (make-array 0 :element-type '(signed-byte 32)) + :type (simple-array (signed-byte 32) (*))) + (roots (make-array 0 :element-type '(signed-byte 32)) + :type (simple-array (signed-byte 32) (*))) + (tree-leaf-offsets (make-array 0 :element-type '(unsigned-byte 32)) + :type (simple-array (unsigned-byte 32) (*)))) + +(defun check-packable (forest) + "Signal PACKED-BUILD-ERROR unless FOREST can be packed." + (let ((dtrees (forest-dtree-list forest)) + (dim (forest-datum-dim forest))) + (unless (= (length dtrees) (forest-n-tree forest)) + (error 'packed-build-error + :detail (format nil "~D trees in the list, FOREST-N-TREE says ~D" + (length dtrees) (forest-n-tree forest)))) + (labels ((walk (node depth) + (cond + ((null node) + (error 'packed-build-error :detail "a nil node")) + ((node-test-attribute node) + (let ((f (node-test-attribute node))) + (unless (and (typep f 'fixnum) (<= 0 f) (< f dim)) + (error 'packed-build-error + :detail (format nil "feature ~S outside [0,~D) at depth ~D" + f dim depth)))) + (unless (typep (node-test-threshold node) 'single-float) + (error 'packed-build-error + :detail (format nil "threshold ~S is not a single-float" + (node-test-threshold node)))) + (unless (and (node-left-node node) (node-right-node node)) + (error 'packed-build-error + :detail "an internal node with only one child")) + (walk (node-left-node node) (1+ depth)) + (walk (node-right-node node) (1+ depth))) + (t + ;; NIL is not the only way to have nothing to count. A split whose + ;; sampled attribute is constant over the node's rows gets + ;; threshold = min = max from MAKE-RANDOM-TEST, and with the >= + ;; convention every row goes left -- leaving the right child a leaf whose + ;; SAMPLE-INDICES is a real but zero-length array. CLASS-DISTRIBUTION + ;; divides by a zero sum either way and returns a uniform distribution + ;; without signalling. Measured: 3 of letter's 26936 leaves at + ;; :max-depth 20. + (let ((indices (node-sample-indices node))) + (unless (and indices (plusp (length indices))) + (error 'packed-build-error + :detail (format nil "a leaf at depth ~D has ~:[no sample ~ +indices at all~;an empty sample-indices array~], so its class distribution would come out ~ +uniform without any error being signalled -- if the forest has been pruned, rebuild it ~ +with :remove-sample-indices? nil (issue #14)" depth indices)))))))) + (dolist (dtree dtrees) + (walk (dtree-root dtree) 0))))) + +(defun count-tree-nodes (dtree) + "Return (values n-internal n-leaf) for DTREE." + (let ((internal 0) (leaves 0)) + (labels ((walk (node) + (if (node-test-attribute node) + (progn (incf internal) + (walk (node-left-node node)) + (walk (node-right-node node))) + (incf leaves)))) + (walk (dtree-root dtree))) + (values internal leaves))) + +(defun check-forward-numbering (n-internal left right) + "Signal PACKED-BUILD-ERROR unless every internal child index exceeds its parent's. + +EMIT produces this by construction. Checking it here costs one pass and turns a later change +to the numbering -- breadth-first, or anything cache-conscious -- into an error at build time +with this message, rather than into files PACKED-LOAD starts rejecting for no visible reason." + (dotimes (i n-internal) + (dolist (side (list (cons "left" (aref left i)) (cons "right" (aref right i)))) + (let ((child (cdr side))) + (when (and (not (minusp child)) (<= child i)) + (error 'packed-build-error + :detail (format nil "~A[~D] = ~D does not exceed its parent index; the ~ +node numbering is no longer depth-first, and PACKED-LOAD's termination check assumes it is" + (car side) i child))))))) + +(defun build-packed-topology (forest) + "Flatten FOREST's structure into arrays. + +Two passes: count and validate, then fill exact-size arrays. Validating first matters +because the traversal runs at (safety 0)." + (check-packable forest) + (let* ((dtrees (forest-dtree-list forest)) + (n-tree (length dtrees)) + (n-internal 0) + (n-leaf 0)) + (dolist (dtree dtrees) + (multiple-value-bind (i l) (count-tree-nodes dtree) + (incf n-internal i) + (incf n-leaf l))) + (unless (< n-internal (expt 2 31)) + (error 'packed-build-error + :detail (format nil "~D internal nodes exceeds a (signed-byte 32) index" + n-internal))) + (unless (< n-leaf (expt 2 31)) + (error 'packed-build-error + :detail (format nil "~D leaves exceeds a (signed-byte 32) index" n-leaf))) + (let ((feature (make-array n-internal :element-type '(unsigned-byte 32))) + (threshold (make-array n-internal :element-type 'single-float)) + (left (make-array n-internal :element-type '(signed-byte 32))) + (right (make-array n-internal :element-type '(signed-byte 32))) + (roots (make-array n-tree :element-type '(signed-byte 32))) + (offsets (make-array n-tree :element-type '(unsigned-byte 32))) + (node-cursor 0) + (leaf-base 0)) + (loop for dtree in dtrees + for tree from 0 + do (setf (aref offsets tree) leaf-base) + (let ((local-leaf 0)) + (labels ((emit (node) + (cond + ((node-test-attribute node) + (let ((self node-cursor)) + (incf node-cursor) + (setf (aref feature self) (node-test-attribute node) + (aref threshold self) (node-test-threshold node)) + (setf (aref left self) (emit (node-left-node node))) + (setf (aref right self) (emit (node-right-node node))) + self)) + (t + (prog1 (lognot (+ leaf-base local-leaf)) + (incf local-leaf)))))) + (setf (aref roots tree) (emit (dtree-root dtree)))) + (incf leaf-base local-leaf))) + (check-forward-numbering n-internal left right) + (%make-packed-topology + :n-tree n-tree :n-internal n-internal :n-leaf n-leaf + :datum-dim (forest-datum-dim forest) + :feature feature :threshold threshold :left left :right right + :roots roots :tree-leaf-offsets offsets)))) + +(declaim (inline check-datamatrix-width)) +(defun check-datamatrix-width (topology datamatrix) + "Signal an error unless DATAMATRIX is at least as wide as the forest was trained on. + +PACKED-LEAF reads DATAMATRIX[datum, FEATURE[node]] at (safety 0), and FEATURE is bounded only +by the training width. Hand the walk a narrower matrix -- the easy mistake, predicting with a +different dataset than you trained on -- and it reads past the row rather than complaining. +This is one comparison against several hundred tree walks per call." + (declare (optimize (speed 3) (safety 0)) + (type packed-topology topology) + (type (simple-array single-float (* *)) datamatrix)) + (let ((dim (packed-topology-datum-dim topology)) + (width (array-dimension datamatrix 1))) + (declare (type fixnum dim width)) + (when (< width dim) + (error "this model was packed from a forest trained on ~D columns, but the datamatrix ~ +has ~D" dim width)))) + +(declaim (inline packed-leaf)) +(defun packed-leaf (topology datamatrix datum-index root) + "The leaf number a datum reaches from ROOT. This is the Global Refinement index." + (declare (optimize (speed 3) (safety 0)) + (type packed-topology topology) + (type (simple-array single-float (* *)) datamatrix) + (type fixnum datum-index) + (type (signed-byte 32) root)) + (let ((feature (packed-topology-feature topology)) + (threshold (packed-topology-threshold topology)) + (left (packed-topology-left topology)) + (right (packed-topology-right topology)) + (node root)) + (declare (type (simple-array (unsigned-byte 32) (*)) feature) + (type (simple-array single-float (*)) threshold) + (type (simple-array (signed-byte 32) (*)) left right) + (type (signed-byte 32) node)) + (loop while (>= node 0) + do (setf node (if (>= (aref datamatrix datum-index (aref feature node)) + (aref threshold node)) + (aref left node) + (aref right node)))) + (lognot node))) + +(defun packed-leaf-indices (topology datamatrix datum-index out) + "Fill OUT with the leaf each tree sends a datum to. OUT is (simple-array fixnum (n-tree))." + (declare (optimize (speed 3) (safety 0)) + (type packed-topology topology) + (type (simple-array single-float (* *)) datamatrix) + (type (simple-array fixnum (*)) out) + (type fixnum datum-index)) + (let ((roots (packed-topology-roots topology)) + (n-tree (packed-topology-n-tree topology))) + (declare (type (simple-array (signed-byte 32) (*)) roots) + (type fixnum n-tree)) + (check-datamatrix-width topology datamatrix) + (dotimes (tree n-tree out) + (setf (aref out tree) + (packed-leaf topology datamatrix datum-index (aref roots tree)))))) diff --git a/src/random-forest.lisp b/src/random-forest.lisp index edb61f5..82fa024 100644 --- a/src/random-forest.lisp +++ b/src/random-forest.lisp @@ -1,6 +1,41 @@ (defpackage :cl-random-forest/src/random-forest (:use #:cl #:cl-random-forest/src/utils) + ;; CLOL, CLOL.VECTOR and ALEXANDRIA are all referenced qualified below (CLOL:MAKE-ONE-VS-REST + ;; and so on) but were never named here, so ASDF's package-inferred-system never loaded them + ;; on a path that does not go through the top-level cl-random-forest system -- exactly the + ;; same defect fixed in src/utils.lisp's DEFPACKAGE, just one file further down the chain that + ;; cl-random-forest/src/packed pulls in. CLOL and CLOL.VECTOR need the extra + ;; ASDF:REGISTER-SYSTEM-PACKAGES mapping in cl-random-forest.asd for the same reason SVMFORMAT + ;; does: both are nicknames of packages provided by CL-ONLINE-LEARNING, not systems of their + ;; own. Importing these symbols does not change how the code below reads or behaves -- they + ;; stay qualified at every call site -- it only tells ASDF what to load first. + (:import-from #:clol + #:copy-learner + #:make-one-vs-rest + #:make-sparse-arow + #:make-sparse-rls + #:one-vs-rest-predict + #:sparse-arow-predict + #:sparse-arow-update + #:sparse-rls-predict + #:sparse-rls-update + #:one-vs-rest-input-dimension + #:one-vs-rest-learner-activate + #:one-vs-rest-learner-bias + #:one-vs-rest-learners-vector + #:one-vs-rest-learner-update + #:one-vs-rest-learner-weight + #:one-vs-rest-n-class + #:sparse-arow-input-dimension + #:sparse-arow-weight) + (:import-from #:clol.vector + #:make-sparse-vector + #:sparse-vector-index-vector + #:sparse-vector-value-vector) + (:import-from #:alexandria + #:iota + #:positive-integer) (:export #:make-dtree #:predict-dtree #:test-dtree diff --git a/src/utils.lisp b/src/utils.lisp index 460120d..5a69f20 100644 --- a/src/utils.lisp +++ b/src/utils.lisp @@ -1,6 +1,22 @@ (defpackage :cl-random-forest/src/utils (:use :cl) (:nicknames :cl-random-forest.utils :clrf.utils) + ;; LPARALLEL, ALEXANDRIA and SVMFORMAT are all used below (the parallelisation macros, + ;; PUSH-NTIMES's own macro-defining code, and the libsvm-format readers, respectively), + ;; but none of them were named here. Naming them is what tells ASDF's package-inferred- + ;; system to load them: without this, loading any system that reaches this file without + ;; going through the top-level cl-random-forest system -- cl-random-forest/src/packed, + ;; for instance -- fails on a cold cache with "Package LPARALLEL does not exist", then, + ;; once that is fixed, "Package SVMFORMAT does not exist". + ;; + ;; SVMFORMAT needs one more thing beyond this clause: it is a nickname of the + ;; CL-LIBSVM-FORMAT package, and package-inferred-system derives a dependency's system + ;; name by downcasing the package name, which would look for a nonexistent "svmformat" + ;; system. cl-random-forest.asd registers the real mapping with + ;; ASDF:REGISTER-SYSTEM-PACKAGES before this file is ever read. + (:import-from #:lparallel #:*kernel*) + (:import-from #:alexandria #:with-gensyms #:once-only) + (:import-from #:svmformat #:parse-file) (:export #:random-uniform #:random-normal #:dotimes/pdotimes diff --git a/t/packed.lisp b/t/packed.lisp new file mode 100644 index 0000000..1adc0dd --- /dev/null +++ b/t/packed.lisp @@ -0,0 +1,411 @@ +(in-package :cl-user) + +(defpackage cl-random-forest-test/packed + (:use :cl :rove :cl-random-forest :cl-random-forest/src/packed + :cl-random-forest-test/fixture) + (:import-from #:cl-random-forest/src/random-forest + #:find-leaf + #:dtree-root + #:node-leaf-index + #:node-test-attribute + #:node-left-node + #:node-sample-indices)) +(in-package :cl-random-forest-test/packed) + +(defun synthetic-forest (&key (remove-sample-indices? nil) (n-tree 30)) + "A small forest on the fixture's deterministic classification set." + (multiple-value-bind (datamatrix target) (synthetic-classification-train) + (make-forest +synthetic-n-class+ datamatrix target + :n-tree n-tree :bagging-ratio 0.3 :max-depth 7 :n-trial 10 + :remove-sample-indices? remove-sample-indices?))) + +(deftest packed-topology-leaf-numbers-are-the-refine-index + ;; The leaf number a topology reports must be the index Global Refinement uses, so a + ;; refine dataset can be built straight from it. The spec makes this a guarantee rather + ;; than a coincidence of two traversal orders agreeing. + (with-serial-kernel + (multiple-value-bind (datamatrix target) (synthetic-classification-train) + (declare (ignore target)) + (let* ((forest (synthetic-forest)) + (topology (build-packed-topology forest)) + (offsets (forest-index-offset forest)) + (out (make-array (forest-n-tree forest) :element-type 'fixnum)) + (checked 0) + (mismatch 0)) + (dotimes (i 200) + (packed-leaf-indices topology datamatrix i out) + (loop for dtree in (forest-dtree-list forest) + for tree from 0 + do (incf checked) + (unless (= (aref out tree) + (+ (node-leaf-index + (find-leaf (dtree-root dtree) datamatrix i)) + (aref offsets tree))) + (incf mismatch)))) + (ok (plusp checked) (format nil "checked ~D (datum, tree) pairs" checked)) + (ok (zerop mismatch) + (format nil "~D leaf numbers differ from the refine index" mismatch)))))) + +(deftest packed-topology-refuses-a-forest-with-stranded-leaves + ;; CLASS-DISTRIBUTION returns a *uniform* distribution for a leaf with no sample indices + ;; rather than signalling, so a builder that read one would freeze that into the model + ;; with nothing to show for it. + ;; + ;; Pruning is what creates such a leaf, not forest construction. SET-BEST-CHILDREN! nils + ;; the indices of the node it *splits*, and DELETE-CHILDREN! later turns that node back + ;; into a leaf without restoring them (issue #14). A leaf is never split, so it keeps its + ;; indices whatever :remove-sample-indices? says -- which is why this test prunes. + (with-serial-kernel + (multiple-value-bind (datamatrix target) (synthetic-classification-train) + (let* ((forest (make-forest +synthetic-n-class+ datamatrix target + :n-tree 30 :bagging-ratio 0.3 :max-depth 7 :n-trial 10 + :remove-sample-indices? t)) + (refine-dataset (make-refine-dataset forest datamatrix)) + (learner (make-refine-learner forest))) + (ok (packed-topology-p (build-packed-topology forest)) + "before pruning, every leaf still has its indices and the forest packs") + (dotimes (epoch 3) + (train-refine-learner learner refine-dataset target)) + (pruning! forest learner 0.3) + (ok (handler-case (progn (build-packed-topology forest) nil) + (packed-build-error () t)) + "after pruning, build-packed-topology signals on the stranded leaves"))))) + +(deftest packed-topology-refuses-a-leaf-with-no-samples-to-count + ;; A leaf can have nothing to count without its SAMPLE-INDICES being NIL. When a split's + ;; sampled attribute is constant over the node's rows, MAKE-RANDOM-TEST produces + ;; threshold = min = max and every row goes left, leaving the right child a leaf holding + ;; a real but zero-length array. CLASS-DISTRIBUTION divides by a zero sum in both cases + ;; and returns a uniform distribution rather than signalling, so both must be rejected. + ;; It is rare but not hypothetical: 3 of letter's 26936 leaves at :max-depth 20. + (with-serial-kernel + (multiple-value-bind (datamatrix target) (synthetic-classification-train) + (declare (ignore target)) + (declare (ignorable datamatrix)) + (let ((forest (synthetic-forest))) + (ok (packed-topology-p (build-packed-topology forest)) + "the forest packs while every leaf has samples") + ;; Empty one leaf's indices, which is what such a split leaves behind. + (labels ((leftmost-leaf (node) + (if (node-test-attribute node) + (leftmost-leaf (node-left-node node)) + node))) + (setf (node-sample-indices + (leftmost-leaf (dtree-root (first (forest-dtree-list forest))))) + (make-array 0 :element-type 'fixnum))) + (ok (handler-case (progn (build-packed-topology forest) nil) + (packed-build-error () t)) + "a leaf with an empty sample-indices array is rejected"))))) + +(deftest packed-topology-counts-add-up + (with-serial-kernel + (let* ((forest (synthetic-forest)) + (topology (build-packed-topology forest))) + (ok (= (packed-topology-n-tree topology) (forest-n-tree forest)) + "one root per tree") + ;; A binary tree with L leaves has L-1 internal nodes, so a forest of T trees has + ;; n-leaf - T of them. + (ok (= (packed-topology-n-internal topology) + (- (packed-topology-n-leaf topology) (packed-topology-n-tree topology))) + (format nil "~D internal, ~D leaves, ~D trees" + (packed-topology-n-internal topology) + (packed-topology-n-leaf topology) + (packed-topology-n-tree topology)))))) + +(deftest packed-classifier-agrees-bit-for-bit + ;; Agreement on the class alone is too weak -- two different distributions can share an + ;; argmax. Compare the whole distribution, and require exact equality: the same floats + ;; are summed in the same tree and class order, so anything else is a defect. + (with-serial-kernel + (multiple-value-bind (datamatrix target) (synthetic-classification-test) + (declare (ignore target)) + (let ((forest (synthetic-forest))) + (dolist (payload '(:dense :csr)) + (let ((classifier (build-packed-classifier forest :payload payload))) + (multiple-value-bind (class-bad dist-bad worst) + (packed-verify classifier forest datamatrix) + (ok (zerop class-bad) + (format nil "~A: ~D classes differ" payload class-bad)) + (ok (zerop dist-bad) + (format nil "~A: ~D distributions differ" payload dist-bad)) + (ok (zerop worst) + (format nil "~A: worst absolute difference ~,10F" payload worst))))))))) + +(deftest packed-classifier-auto-picks-on-the-dense-row-size + ;; CSR once a dense row exceeds a cache line, dense below. Not the compression ratio: + ;; a 10-class model compresses 2.4x on that measure and runs 0.82x as fast. + ;; + ;; The rule is unit-tested on both sides of the boundary, because no fixture forest has + ;; enough classes to reach the CSR side -- the synthetic set has 4. + (let ((rule #'cl-random-forest/src/packed/classifier::choose-payload)) + (ok (eq :dense (funcall rule 1)) "1 class is dense") + (ok (eq :dense (funcall rule 16)) "16 classes is dense, a row being exactly 64 bytes") + (ok (eq :csr (funcall rule 17)) "17 classes is csr, a row exceeding a cache line") + (ok (eq :csr (funcall rule 26)) "26 classes is csr")) + (with-serial-kernel + (let ((forest (synthetic-forest))) + (ok (eq :dense (packed-classifier-kind (build-packed-classifier forest))) + (format nil "a ~D-class forest selects dense" +synthetic-n-class+))))) + +(deftest packed-classifier-csr-holds-only-the-non-zero-entries + (with-serial-kernel + (let* ((forest (synthetic-forest)) + (topology (build-packed-topology forest)) + (dense (build-packed-classifier forest :payload :dense)) + (csr (build-packed-classifier forest :payload :csr)) + (n-leaf (packed-topology-n-leaf topology)) + (n-class (packed-classifier-n-class dense)) + (expected 0)) + (dotimes (row n-leaf) + (dotimes (k n-class) + (unless (zerop (aref (packed-classifier-table dense) row k)) + (incf expected)))) + (ok (= expected (length (packed-classifier-probability csr))) + (format nil "~D non-zero entries, CSR stores ~D" + expected (length (packed-classifier-probability csr))))))) + +(deftest packed-batch-agrees-with-per-datum + ;; Tree-major batching walks one tree over a whole tile before moving to the next, which + ;; changes the order arrays are touched but not the arithmetic: each datum's accumulator + ;; receives the same values in the same tree order. + (with-serial-kernel + (multiple-value-bind (datamatrix target) (synthetic-classification-test) + (declare (ignore target)) + (let* ((forest (synthetic-forest)) + (n-rows (array-dimension datamatrix 0))) + (dolist (payload '(:dense :csr)) + (let* ((classifier (build-packed-classifier forest :payload payload)) + (acc (make-packed-accumulator classifier)) + (bad 0)) + (dolist (tile (list 1 7 64 n-rows)) + (let ((accs (make-packed-accumulators classifier tile)) + (out (make-array tile :element-type 'fixnum)) + (i 0)) + (loop while (< i n-rows) + do (let ((end (min n-rows (+ i tile)))) + (packed-predict-batch classifier datamatrix i end accs out) + (loop for j from i below end + do (unless (= (aref out (- j i)) + (packed-predict classifier datamatrix j acc)) + (incf bad))) + (setf i end))))) + (ok (zerop bad) + (format nil "~A: ~D batch answers differ from per-datum" payload bad)))))))) + +(deftest packed-round-trips-through-a-file + (with-serial-kernel + (multiple-value-bind (datamatrix target) (synthetic-classification-test) + (declare (ignore target)) + (let ((forest (synthetic-forest))) + (dolist (payload '(:dense :csr)) + (uiop:with-temporary-file (:pathname path :type "packed") + (let* ((original (build-packed-classifier forest :payload payload)) + (acc (make-packed-accumulator original))) + (packed-save original path) + (let* ((restored (packed-load path)) + (restored-acc (make-packed-accumulator restored)) + (bad 0)) + (ok (eq (packed-classifier-kind restored) payload) + (format nil "~A survives the round trip" payload)) + (dotimes (i (array-dimension datamatrix 0)) + (unless (= (packed-predict original datamatrix i acc) + (packed-predict restored datamatrix i restored-acc)) + (incf bad))) + (ok (zerop bad) + (format nil "~A: ~D predictions differ after reload" payload bad)))))))))) + +(deftest packed-load-rejects-a-file-it-cannot-trust + (with-serial-kernel + (let ((forest (synthetic-forest))) + (uiop:with-temporary-file (:pathname path :type "packed") + (packed-save (build-packed-classifier forest) path) + ;; Corrupt the magic. + (with-open-file (s path :direction :io :if-exists :overwrite + :element-type '(unsigned-byte 8)) + (file-position s 0) + (write-byte 0 s)) + (ok (handler-case (progn (packed-load path) nil) + (packed-load-error () t)) + "a bad magic number is rejected")) + (uiop:with-temporary-file (:pathname path :type "packed") + (packed-save (build-packed-classifier forest) path) + ;; Corrupt the version, which lives at byte offset 8. + (with-open-file (s path :direction :io :if-exists :overwrite + :element-type '(unsigned-byte 8)) + (file-position s 8) + (write-byte 99 s)) + (ok (handler-case (progn (packed-load path) nil) + (packed-load-error () t)) + "an unknown version is rejected")) + (uiop:with-temporary-file (:pathname path :type "packed") + (packed-save (build-packed-classifier forest) path) + ;; Corrupt the byte-order probe, at offset 12. + (with-open-file (s path :direction :io :if-exists :overwrite + :element-type '(unsigned-byte 8)) + (file-position s 12) + (write-byte 99 s)) + (ok (handler-case (progn (packed-load path) nil) + (packed-load-error () t)) + "a byte-order mismatch is rejected")) + (uiop:with-temporary-file (:pathname path :type "packed") + (packed-save (build-packed-classifier forest) path) + ;; Truncate to half length, keeping the surviving bytes intact. Overwriting with + ;; zeros instead (as an earlier version of this test did) zeros the magic along + ;; with everything else, so PACKED-LOAD would reject the file on the magic check + ;; the first case above already covers, and the short-read paths this case exists + ;; to exercise would never run. + (let* ((bytes (with-open-file (s path :element-type '(unsigned-byte 8)) + (file-length s))) + (buffer (make-array bytes :element-type '(unsigned-byte 8)))) + (with-open-file (s path :element-type '(unsigned-byte 8)) + (read-sequence buffer s)) + (with-open-file (s path :direction :output :if-exists :supersede + :element-type '(unsigned-byte 8)) + (write-sequence (subseq buffer 0 (floor bytes 2)) s))) + (ok (handler-case (progn (packed-load path) nil) + (packed-load-error () t)) + "a truncated file is rejected")) + (uiop:with-temporary-file (:pathname path :type "packed") + ;; Corrupt N-CLASS, the seventh of the header's eight u32 fields, at byte offset + ;; 8 + 4*6 = 32. Decrementing it below the real class count is exactly what the + ;; reviewer demonstrated slipping through unchecked: the header positivity check + ;; still passes, and the topology is untouched, but a CSR file's CLASS array now + ;; holds ids that no longer fit under the corrupted count, which + ;; VALIDATE-CSR-ARRAYS must catch. Forcing :CSR here (rather than the :AUTO the + ;; other cases use) is what puts a CLASS array in the file for the corruption to + ;; land in -- the fixture's 4-class forest packs dense under :AUTO. + (packed-save (build-packed-classifier forest :payload :csr) path) + (with-open-file (s path :direction :io :if-exists :overwrite + :element-type '(unsigned-byte 8)) + (file-position s 32) + (write-byte 2 s)) + (ok (handler-case (progn (packed-load path) nil) + (packed-load-error () t)) + "a corrupted n-class is rejected")) + (flet ((clobber (payload offset bytes) + "Save FOREST, overwrite BYTES at OFFSET, and return T if the load is refused." + (uiop:with-temporary-file (:pathname path :type "packed") + (packed-save (build-packed-classifier forest :payload payload) path) + (with-open-file (s path :direction :io :if-exists :overwrite + :element-type '(unsigned-byte 8)) + (file-position s offset) + (dolist (b bytes) (write-byte b s))) + (handler-case (progn (packed-load path) nil) + (packed-load-error () t)))) + (u32-at (offset) + "Read the u32 the saved header holds at OFFSET, to avoid hardcoding counts." + (uiop:with-temporary-file (:pathname path :type "packed") + (packed-save (build-packed-classifier forest) path) + (with-open-file (s path :element-type '(unsigned-byte 8)) + (file-position s offset) + (+ (read-byte s) (ash (read-byte s) 8) + (ash (read-byte s) 16) (ash (read-byte s) 24)))))) + ;; The topology arrays begin right after the header's eight u32 fields, at byte + ;; offset 8 + 4*8 = 40. FEATURE comes first and THRESHOLD second, so LEFT starts at + ;; 40 + 4*n-internal*2. N-INTERNAL is read back from the header (offset 24) rather + ;; than hardcoded, since the exact node count is an artifact of training the fixture + ;; forest, not something this test pins. + (let* ((n-internal (u32-at 24)) + (feature-start 40) + (left-start (+ 40 (* 4 n-internal 2)))) + (ok (clobber :auto left-start '(#xff #xff #xff #x7f)) + "a corrupted left entry is rejected") + ;; A node whose left child is itself. Every bound holds -- 0 is a perfectly good + ;; internal-node index -- so only the forward-motion rule catches it. Without that + ;; rule PACKED-LEAF spins in this node forever at (safety 0), which is why the + ;; check exists: a hang is not something a caller can recover from. + (ok (clobber :auto left-start '(0 0 0 0)) + "a left entry pointing at its own node is rejected") + ;; A feature index past the datamatrix's columns. Nothing about the tree structure + ;; is wrong here; the walk would simply read off the end of a row. + (ok (clobber :auto feature-start '(#xff #xff #xff #xff)) + "a feature index outside the trained datamatrix width is rejected") + ;; An inflated N-INTERNAL. READ-ARRAY has to size its buffers from this count + ;; before it can notice the file is too short, so without the remaining-length + ;; check this asks for 16 GB and dies of heap exhaustion rather than signalling. + ;; The assertion is specifically that it is a PACKED-LOAD-ERROR. + (ok (clobber :auto 24 '(#xff #xff #xff #xff)) + "a count larger than the file is rejected before anything is allocated")))))) + +(deftest packed-refuses-a-datamatrix-narrower-than-the-forest + ;; FEATURE is bounded by the training width, and the walk reads the datamatrix at + ;; (safety 0), so a narrower matrix is an out-of-bounds read rather than an error. This is + ;; reachable with a perfectly good in-memory model -- no file involved -- by predicting + ;; with a different dataset than the forest was trained on. + (with-serial-kernel + (let* ((forest (synthetic-forest)) + (classifier (build-packed-classifier forest)) + (dim (packed-topology-datum-dim (packed-classifier-topology classifier))) + (acc (make-packed-accumulator classifier)) + (narrow (make-array (list 1 (1- dim)) :element-type 'single-float + :initial-element 0.0))) + (ok (plusp dim) "the packed topology remembers the training width") + (ok (handler-case (progn (packed-predict classifier narrow 0 acc) nil) + (error () t)) + "predicting from a too-narrow datamatrix signals instead of reading past the row")))) + +(deftest packed-float-bits-round-trip + ;; The fast paths exist per implementation; the portable fallback is the reference. They + ;; must agree, or a model saved on one implementation would not load on another. + ;; + ;; Both directions are cross-checked, and the list includes true denormals. An earlier + ;; version of this test used only LEAST-POSITIVE-NORMALIZED-SINGLE-FLOAT and compared + ;; encoding alone, and missed that the portable encoder wrote a bogus exponent field for + ;; every denormal -- 1.4012985e-45 came out as 00800001 where IEEE-754 says 00000001. + (let ((values (list 0.0 -0.0 1.0 -1.0 0.5 -0.5 3.14159 1.0e-8 1.0e8 + most-positive-single-float + least-positive-normalized-single-float + (/ least-positive-normalized-single-float 2.0) + least-positive-single-float + (- least-positive-single-float) + (* 12345.0 least-positive-single-float)))) + (dolist (v values) + (ok (= v (cl-random-forest/src/packed::bits-to-single-float + (cl-random-forest/src/packed::single-float-to-bits v))) + (format nil "~S survives the bit round trip" v)) + (ok (= (cl-random-forest/src/packed::single-float-to-bits v) + (cl-random-forest/src/packed::%portable-single-float-to-bits v)) + (format nil "~S: the fast encoder agrees with the portable one" v)) + (ok (= v (cl-random-forest/src/packed::%portable-bits-to-single-float + (cl-random-forest/src/packed::single-float-to-bits v))) + (format nil "~S: the portable decoder agrees with the fast encoder" v))))) + +(deftest packed-float-bits-do-not-depend-on-the-decode-convention + ;; INTEGER-DECODE-FLOAT owes its caller only a (significand, exponent) pair whose product + ;; is the magnitude. Which member of that equivalence class it picks is up to the + ;; implementation, and for denormals the two we run on disagree: SBCL leaves the + ;; significand unnormalised (LEAST-POSITIVE-SINGLE-FLOAT is 1 x 2^-149) while CCL + ;; normalises it (2^23 x 2^-172). + ;; + ;; The portable encoder used to decide "is this a denormal?" by asking whether the + ;; significand was below 2^23 -- reading the host's convention rather than the format. + ;; That was right on SBCL and wrong on CCL, where 5.877472e-39 encoded as 00000000: a + ;; denormal silently becoming zero. Only CI caught it, because on SBCL the host never + ;; produces the pair that breaks it. + ;; + ;; So this feeds the pairs in directly. It fails on any implementation if the convention + ;; ever leaks back into the arithmetic. + (labels ((normalised (s e) + "The CCL convention: significand shifted up into [2^23, 2^24)." + (loop while (< s (ash 1 23)) do (setf s (ash s 1)) (decf e)) + (values s e)) + (minimal (s e) + "The other extreme: significand reduced until odd." + (loop while (and (plusp s) (evenp s)) do (setf s (ash s -1)) (incf e)) + (values s e))) + (dolist (v (list least-positive-single-float + (- least-positive-single-float) + (* 12345.0 least-positive-single-float) + (/ least-positive-normalized-single-float 2.0) + least-positive-normalized-single-float + 1.0 -0.5 3.14159 most-positive-single-float)) + (multiple-value-bind (s e sign) (integer-decode-float v) + (let ((expected (cl-random-forest/src/packed::single-float-to-bits v))) + (dolist (convention (list (cons "as decoded" (list s e)) + (cons "normalised" (multiple-value-list (normalised s e))) + (cons "minimal" (multiple-value-list (minimal s e))) + (cons "over-scaled" (list (ash s 5) (- e 5))))) + (ok (= expected + (cl-random-forest/src/packed::%float-parts-to-bits + (first (cdr convention)) (second (cdr convention)) sign)) + (format nil "~S: ~A significand encodes the same" v (car convention)))))))))