Repository: chrovis/cljam
Pull request #330 (2024-12-24 02:06:35)

rANS (rans4x8) encoding support for CRAM block compression
This PR implements the rANS (rans4x8) encoder to allow the CRAM writer to use it for block compression.

With this PR, the default compressor for the following data series will be updated to rans4x8:

| Data Series | Default Compressor |
| :---: | :---: |
| `BF` | rans4x8 (Order 1) |
| `CF` | rans4x8 (Order 1) |
| `RI` | rans4x8 (Order 0) |
| `RL` | rans4x8 (Order 1) |
| `AP` | rans4x8 (Order 0) |
| `RG` | rans4x8 (Order 1) |
| `NS` | rans4x8 (Order 1) |
| `TS` | rans4x8 (Order 1) |
| `BA` | rans4x8 (Order 1) |
| `QS` | rans4x8 (Order 1) |

You can override the compressor for specific data series or tags by specifying `:r4x8-o0` (Order 0) or `:r4x8-o1` (Order 1) in the `:ds-compressor-overrides` / `:tag-compressor-overrides` options. 

For example, to set the compressor of the `FN` data series to rANS (Order 0) and the `TL` data series to rANS (Order 1), you can configure the writer as follows:

```clojure
(require '[cljam.io.cram :as cram])

(with-open [w (cram/writer "path/to/cram/file"
                           { ...
                            :ds-compressor-overrides {:FN :r4x8-o0
                                                      :TL :r4x8-o1}
                             ... })]
  ... )


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/chrovis/cljam/pull/330?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis) Report
Attention: Patch coverage is `84.64164%` with `45 lines` in your changes missing coverage. Please review.
> Project coverage is 89.81%. Comparing base [(`18f8cca`)](https://app.codecov.io/gh/chrovis/cljam/commit/18f8cca90f5f210f0239dfce2134fbc02e8c044c?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis) to head [(`70fea79`)](https://app.codecov.io/gh/chrovis/cljam/commit/70fea79899f2fa69065dbab19b2e3489d0875825?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis).

| [Files with missing lines](https://app.codecov.io/gh/chrovis/cljam/pull/330?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis) | Patch % | Lines |
|---|---|---|
| [src/cljam/io/cram/codecs/rans4x8.clj](https://app.codecov.io/gh/chrovis/cljam/pull/330?src=pr&el=tree&filepath=src%2Fcljam%2Fio%2Fcram%2Fcodecs%2Frans4x8.clj&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis#diff-c3JjL2NsamFtL2lvL2NyYW0vY29kZWNzL3JhbnM0eDguY2xq) | 84.55% | [4 Missing and 38 partials :warning: ](https://app.codecov.io/gh/chrovis/cljam/pull/330?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis) |
| [src/cljam/io/cram/encode/compressor.clj](https://app.codecov.io/gh/chrovis/cljam/pull/330?src=pr&el=tree&filepath=src%2Fcljam%2Fio%2Fcram%2Fencode%2Fcompressor.clj&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis#diff-c3JjL2NsamFtL2lvL2NyYW0vZW5jb2RlL2NvbXByZXNzb3IuY2xq) | 80.00% | [2 Missing :warning: ](https://app.codecov.io/gh/chrovis/cljam/pull/330?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis) |
| [src/cljam/io/cram/encode/structure.clj](https://app.codecov.io/gh/chrovis/cljam/pull/330?src=pr&el=tree&filepath=src%2Fcljam%2Fio%2Fcram%2Fencode%2Fstructure.clj&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis#diff-c3JjL2NsamFtL2lvL2NyYW0vZW5jb2RlL3N0cnVjdHVyZS5jbGo=) | 0.00% | [0 Missing and 1 partial :warning: ](https://app.codecov.io/gh/chrovis/cljam/pull/330?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master     #330      +/-   ##
==========================================
- Coverage   89.97%   89.81%   -0.17%     
==========================================
  Files         104      104              
  Lines        9341     9623     +282     
  Branches      490      528      +38     
==========================================
+ Hits         8405     8643     +238     
- Misses        446      452       +6     
- Partials      490      528      +38     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/chrovis/cljam/pull/330?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=chrovis).

athos: Thank you for the comment!

I hadn’t noticed the overflow issue with the `byte` fn because I didn’t disable `*unchecked-math*` during development. Commit 46b2f739920b987beb9050304bba829873e90a57 addressed this issue for the `cljam.io.cram.codecs.rans4x8` ns (including the decoder), but it seems like many other namespaces also have similar issues, most of which are likely due to `cljam.io.util.byte-buffer` and `cljam.io.util.lsb.*`. So, I would like to handle those in a separate PR, if possible.

<details><summary>Failed tests</summary>

```console
$ lein test
...
Ran 372 tests containing 3830 assertions.
1 failures, 17 errors.
$ lein test | grep 'lein test :only'
lein test :only cljam.algo.cram-indexer-test/cram-without-alignments-test
lein test :only cljam.io.cram-test/writer-with-reference-embedding-test
lein test :only cljam.io.cram-test/writer-test
lein test :only cljam.io.cram-test/writer-index-options-test
lein test :only cljam.io.cram.data-series-test/build-tag-encoders-test
lein test :only cljam.io.cram.data-series-test/build-tag-decoders-test
lein test :only cljam.io.cram.encode.record-test/encode-slice-records-test
lein test :only cljam.io.cram.encode.structure-test/encode-block-test
lein test :only cljam.io.cram.encode.structure-test/encode-slice-header-block-test
lein test :only cljam.io.cram.encode.structure-test/encode-eof-container-test
lein test :only cljam.io.cram.encode.structure-test/encode-cram-header-container-test
lein test :only cljam.io.cram.encode.structure-test/encode-compression-header-block-test
lein test :only cljam.io.cram.itf8-test/encode-itf8-test
lein test :only cljam.io.util-test/writer-predicates-test
lein test :only cljam.io.util.byte-buffer-test/allocate-msb-byte-buffer-test
lein test :only cljam.io.util.byte-buffer-test/allocate-lsb-byte-buffer-test
lein test :only cljam.io.util.byte-buffer-test/allocate-lsb-byte-buffer-test
lein test :only cljam.io.util.lsb.io-stream-test/about-data-output-stream
$
```
</details>
alumi: Thank you for updating! The commit 46b2f73 perfectly fixes the issue! 👍

> many other namespaces also have similar issues, most of which are likely due to cljam.io.util.byte-buffer and cljam.io.util.lsb.*. So, I would like to handle those in a separate PR, if possible.

I really appreciate you taking care of this in a separate PR! ❤️ 
athos: It seems like the CI failed because Leiningen has been removed from the latest `ubuntu-latest` image 🤔 

cf. https://github.com/actions/runner-images/issues/10636
alumi: Ah, that's a bit of a hassle! 😲 How about using [`DeLaGuardo/setup-clojure`](https://github.com/DeLaGuardo/setup-clojure) to install Leiningen?

Patch:
diff --git a/src/cljam/io/cram/codecs/rans4x8.clj b/src/cljam/io/cram/codecs/rans4x8.clj
index ff32a23..7a9846d 100644
--- a/src/cljam/io/cram/codecs/rans4x8.clj
+++ b/src/cljam/io/cram/codecs/rans4x8.clj
@@ -1,7 +1,8 @@
 (ns cljam.io.cram.codecs.rans4x8
-  (:require [cljam.io.util.byte-buffer :as bb]
-            [cljam.io.cram.itf8 :as itf8])
-  (:import [java.util Arrays]))
+  (:require [cljam.io.cram.itf8 :as itf8]
+            [cljam.io.util.byte-buffer :as bb])
+  (:import [java.nio Buffer ByteBuffer]
+           [java.util Arrays]))
 
 (def ^:private byte-array-type (type (byte-array 0)))
 (def ^:private int-array-type (type (int-array 0)))
@@ -49,9 +50,9 @@
         (let [curr (aget cum-freqs i)]
           (if (= start curr)
             (recur (inc i) start)
-            (do (Arrays/fill arr start curr (byte (dec i)))
+            (do (Arrays/fill arr start curr (unchecked-byte (dec i)))
                 (recur (inc i) curr))))
-        (Arrays/fill arr start 4096 (byte 255))))
+        (Arrays/fill arr start 4096 (unchecked-byte 255))))
     arr))
 
 (defn- advance-step ^long [^long c ^long f ^long state]
@@ -79,7 +80,7 @@
             state' (->> state
                         (advance-step (aget cum-freqs sym) (aget freqs sym))
                         (renormalize-state bb))]
-        (aset out i (byte sym))
+        (aset out i (unchecked-byte sym))
         (aset states j state')))
     out))
 
@@ -107,7 +108,7 @@
                           (advance-step (aget cfreqs sym)
                                         (aget ^ints (aget freqs last-sym) sym))
                           (renormalize-state bb))]
-          (aset out (+ i (* j quarter)) (byte sym))
+          (aset out (+ i (* j quarter)) (unchecked-byte sym))
           (aset states j state')
           (aset last-syms j sym))))
     (dotimes [i (- n-out truncated)]
@@ -120,7 +121,7 @@
                         (advance-step (aget cfreq sym)
                                       (aget ^ints (aget freqs last-sym) sym))
                         (renormalize-state bb))]
-        (aset out (+ i truncated) (byte sym))
+        (aset out (+ i truncated) (unchecked-byte sym))
         (aset states 3 state')
         (aset last-syms 3 sym)))
     out))
@@ -135,3 +136,321 @@
     (if (zero? order)
       (decode0 bb n-out)
       (decode1 bb n-out))))
+
+(defn- normalize-frequencies! [^ints freqs ^long total]
+  (let [tr (+ (quot (bit-shift-left 4096 31) total)
+              (quot (bit-shift-left 1 31) total))]
+    (loop [i 0, m 0, M 0, fsum 0]
+      (if (< i 256)
+        (let [f (aget freqs i)]
+          (if (zero? f)
+            (recur (inc i) m M fsum)
+            (let [f' (as-> (unsigned-bit-shift-right (* f tr) 31) f'
+                       (if (zero? f') 1 f'))
+                  fsum' (+ fsum f')]
+              (aset freqs i f')
+              (if (< m f')
+                (recur (inc i) f' i fsum')
+                (recur (inc i) m M fsum')))))
+        (let [f (aget freqs M)
+              fsum' (inc fsum)]
+          (if (< fsum' 4096)
+            (aset freqs M (+ f (- 4096 fsum')))
+            (aset freqs M (- f (- fsum' 4096)))))))))
+
+(defmacro ^:private ainc!
+  [arr i]
+  {:pre [(symbol? arr) (symbol? i)]}
+  `(aset ~arr ~i (inc (aget ~arr ~i))))
+
+(defn- calculate-frequencies0 ^ints [^ByteBuffer bb]
+  (let [total (.remaining bb)
+        freqs (int-array 256)]
+    (dotimes [_ total]
+      (let [b (long (bb/read-ubyte bb))]
+        (ainc! freqs b)))
+    (normalize-frequencies! freqs total)
+    freqs))
+
+(definline ^:private read-ubyte-from [^ByteBuffer bb i]
+  `(bit-and 0xff (.get ~bb (long ~i))))
+
+(defn- calculate-frequencies1 ^"[[I" [^ByteBuffer bb]
+  (let [size (.remaining bb)
+        ^"[[I" freqs (make-array Integer/TYPE 256 256)
+        totals (int-array 256)
+        _ (loop [i 0, prev 0]
+            (when (< i size)
+              (let [b (long (bb/read-ubyte bb))
+                    ^ints fs (aget freqs prev)]
+                (ainc! fs b)
+                (ainc! totals prev)
+                (recur (inc i) b))))
+        q (unsigned-bit-shift-right size 2)
+        ^ints f0 (aget freqs 0)
+        b (read-ubyte-from bb q)
+        _ (ainc! f0 b)
+        b (read-ubyte-from bb (* 2 q))
+        _ (ainc! f0 b)
+        b (read-ubyte-from bb (* 3 q))
+        _ (ainc! f0 b)]
+    (aset totals 0 (+ (aget totals 0) 3))
+    (dotimes [i 256]
+      (let [total (aget totals i)]
+        (when-not (zero? total)
+          (normalize-frequencies! (aget freqs i) total))))
+    freqs))
+
+(defn- next-rle ^long [^long rle ^ints freqs ^long i ^ByteBuffer out]
+  (if (zero? rle)
+    (do (.put out (unchecked-byte i))
+        (if (and (> i 0) (not (zero? (aget freqs (dec i)))))
+          (let [rle' (loop [rle (inc i)]
+                       (if (and (< rle 256)
+                                (not (zero? (aget freqs rle))))
+                         (recur (inc rle))
+                         (- rle (inc i))))]
+            (.put out (unchecked-byte rle'))
+            rle')
+          rle))
+    (dec rle)))
+
+(defn- encode-itf8
+  "Simplified version of ITF8 encoder for up to two bytes"
+  [^ByteBuffer out ^long n]
+  (if (< n 128)
+    (.put out (unchecked-byte n))
+    (do (.put out (unchecked-byte (bit-or 128 (unsigned-bit-shift-right n 8))))
+        (.put out (unchecked-byte (bit-and 0xff n))))))
+
+(defn- write-frequencies0 ^long [^ByteBuffer out ^ints freqs]
+  (let [start (.position out)]
+    (loop [i 0, rle 0]
+      (when (< i 256)
+        (let [f (aget freqs i)]
+          (if (zero? f)
+            (recur (inc i) rle)
+            (let [rle' (next-rle rle freqs i out)]
+              (encode-itf8 out f)
+              (recur (inc i) rle'))))))
+    (.put out (unchecked-byte 0))
+    (- (.position out) start)))
+
+(defn- write-frequencies1 ^long [^ByteBuffer out ^"[[I" freqs]
+  (let [start (.position out)
+        totals (int-array 256)
+        _ (dotimes [i 256]
+            (let [^ints fs (aget freqs i)]
+              (dotimes [j 256]
+                (aset totals i (+ (aget totals i) (aget fs j))))))]
+    (loop [i 0, rle-i 0]
+      (when (< i 256)
+        (if (zero? (aget totals i))
+          (recur (inc i) rle-i)
+          (let [rle-i' (next-rle rle-i totals i out)
+                ^ints fs (aget freqs i)]
+            (loop [j 0, rle-j 0]
+              (when (< j 256)
+                (let [f (aget fs j)]
+                  (if (zero? f)
+                    (recur (inc j) rle-j)
+                    (let [rle-j' (next-rle rle-j fs j out)]
+                      (encode-itf8 out f)
+                      (recur (inc j) rle-j'))))))
+            (.put out (unchecked-byte 0))
+            (recur (inc i) rle-i')))))
+    (.put out (unchecked-byte 0))
+    (- (.position out) start)))
+
+(defprotocol ISymbolState
+  (init! [this start freq])
+  (update! [this ^ByteBuffer bb b]))
+
+(deftype SymbolState
+         [^:unsynchronized-mutable ^long xmax
+          ^:unsynchronized-mutable ^long rcp-freq
+          ^:unsynchronized-mutable ^long bias
+          ^:unsynchronized-mutable ^long cmpl-freq
+          ^:unsynchronized-mutable ^long rcp-shift]
+  ISymbolState
+  (init! [_ start freq]
+    (let [start (long start)
+          freq (long freq)]
+      (set! xmax (* 0x80000 freq))
+      (set! cmpl-freq (- 0x1000 freq))
+      (if (< freq 2)
+        (do (set! rcp-freq (bit-not 0))
+            (set! rcp-shift 0)
+            (set! bias (dec (+ start 0x1000))))
+        (let [shift (long
+                     (loop [shift 0]
+                       (if (< (bit-shift-left 1 shift) freq)
+                         (recur (inc shift))
+                         shift)))]
+          (set! rcp-freq (quot (dec (+ (bit-shift-left 0x80000000 shift) freq)) freq))
+          (set! rcp-shift (dec shift))
+          (set! bias start)))
+      (set! rcp-shift (+ rcp-shift 32))))
+  (update! [_ bb r]
+    (let [x (long
+             (loop [i 2, x (long r)]
+               (if (or (zero? i) (< x xmax))
+                 x
+                 (do (.put ^ByteBuffer bb (unchecked-byte (bit-and 0xff x)))
+                     (recur (dec i) (unsigned-bit-shift-right x 8))))))
+          q (unsigned-bit-shift-right (* x (bit-and 0xffffffff rcp-freq)) rcp-shift)]
+      (+ x bias (* q cmpl-freq)))))
+
+(def ^:private ^:const RANS_BYTE_L 0x800000)
+
+(defn- reverse-buffer! [^ByteBuffer bb]
+  (let [arr (.array bb)
+        offset (.arrayOffset bb)
+        size (.limit bb)]
+    (loop [i offset, j (dec (+ offset size))]
+      (when (< i j)
+        (let [t (aget arr j)]
+          (aset arr j (aget arr i))
+          (aset arr i t)
+          (recur (inc i) (dec j)))))))
+
+(defn- encode-payload0 ^long [^ByteBuffer in ^objects syms ^ByteBuffer out]
+  (let [raw-size (.remaining in)
+        r (bit-and raw-size 3)
+        out' (.slice out)
+        r2 (if (= r 3)
+             (long (update! (aget syms (read-ubyte-from in (- raw-size r -2))) out' RANS_BYTE_L))
+             RANS_BYTE_L)
+        r1 (if (>= r 2)
+             (long (update! (aget syms (read-ubyte-from in (- raw-size r -1))) out' RANS_BYTE_L))
+             RANS_BYTE_L)
+        r0 (if (>= r 1)
+             (long (update! (aget syms (read-ubyte-from in (- raw-size r))) out' RANS_BYTE_L))
+             RANS_BYTE_L)]
+    (loop [i (bit-and raw-size (bit-not 3)), r3 RANS_BYTE_L, r2 r2, r1 r1, r0 r0]
+      (if (> i 0)
+        (let [r3' (long (update! (aget syms (read-ubyte-from in (- i 1))) out' r3))
+              r2' (long (update! (aget syms (read-ubyte-from in (- i 2))) out' r2))
+              r1' (long (update! (aget syms (read-ubyte-from in (- i 3))) out' r1))
+              r0' (long (update! (aget syms (read-ubyte-from in (- i 4))) out' r0))]
+          (recur (- i 4) r3' r2' r1' r0'))
+        (do (.putInt out' r3)
+            (.putInt out' r2)
+            (.putInt out' r1)
+            (.putInt out' r0)
+            (.flip ^Buffer out')
+            (reverse-buffer! out')
+            (.position ^Buffer in (.limit in))
+            (.limit out'))))))
+
+(defmacro ^:private aget2 [syms i j]
+  `(aget ~(with-meta `(aget ~syms ~i) {:tag 'objects}) ~j))
+
+(defn- encode-payload1 ^long [^ByteBuffer in ^objects syms ^ByteBuffer out]
+  (let [raw-size (.remaining in)
+        q (unsigned-bit-shift-right raw-size 2)
+        i0 (- q 2)
+        i1 (- (* 2 q) 2)
+        i2 (- (* 3 q) 2)
+        l0 (if (>= (inc i0) 0) (read-ubyte-from in (inc i0)) 0)
+        l1 (if (>= (inc i1) 0) (read-ubyte-from in (inc i1)) 0)
+        l2 (if (>= (inc i2) 0) (read-ubyte-from in (inc i2)) 0)
+        out' (.slice out)]
+    (loop [i3 (- raw-size 2)
+           l3 (read-ubyte-from in (dec raw-size))
+           r3 RANS_BYTE_L]
+      (if (and (> i3 (- (* 4 q) 2)) (>= i3 0))
+        (let [c3 (read-ubyte-from in i3)
+              r3' (long (update! (aget2 syms c3 l3) out' r3))]
+          (recur (dec i3) c3 r3'))
+        (loop [i0 i0, i1 i1, i2 i2, i3 i3,
+               l0 l0, l1 l1, l2 l2, l3 l3,
+               r0 RANS_BYTE_L, r1 RANS_BYTE_L, r2 RANS_BYTE_L, r3 r3]
+          (if (>= i0 0)
+            (let [c0 (read-ubyte-from in i0)
+                  c1 (read-ubyte-from in i1)
+                  c2 (read-ubyte-from in i2)
+                  c3 (read-ubyte-from in i3)
+                  r3' (long (update! (aget2 syms c3 l3) out' r3))
+                  r2' (long (update! (aget2 syms c2 l2) out' r2))
+                  r1' (long (update! (aget2 syms c1 l1) out' r1))
+                  r0' (long (update! (aget2 syms c0 l0) out' r0))]
+              (recur (dec i0) (dec i1) (dec i2) (dec i3) c0 c1 c2 c3 r0' r1' r2' r3'))
+            (let [r3' (long (update! (aget2 syms 0 l3) out' r3))
+                  r2' (long (update! (aget2 syms 0 l2) out' r2))
+                  r1' (long (update! (aget2 syms 0 l1) out' r1))
+                  r0' (long (update! (aget2 syms 0 l0) out' r0))]
+              (.putInt out' r3')
+              (.putInt out' r2')
+              (.putInt out' r1')
+              (.putInt out' r0')
+              (.flip ^Buffer out')
+              (reverse-buffer! out')
+              (.position ^Buffer in (.limit in))
+              (.limit out'))))))))
+
+(def ^:private ^:const PREFIX_BYTE_LEN (+ 1 4 4))
+
+(defn- allocate-output-buffer [^long raw-size]
+  ;; The size estimation code comes from:
+  ;; - https://github.com/samtools/htscodecs/blob/51794289ac47455209c333182b6768f99a613948/htscodecs/rANS_static.c#L77
+  ;; - https://github.com/samtools/htscodecs/blob/51794289ac47455209c333182b6768f99a613948/htscodecs/rANS_static.c#L410
+  (let [allocated-size (+ (* 1.05 raw-size)
+                          ;; upper bound of frequency table size
+                          (* 257 257 3)
+                          ;; prefix
+                          PREFIX_BYTE_LEN)]
+    (bb/allocate-lsb-byte-buffer allocated-size)))
+
+(defn- encode0 ^long [^Buffer in ^ByteBuffer out]
+  (let [freqs (calculate-frequencies0 in)
+        syms (object-array 256)
+        _ (loop [i 0, total 0]
+            (when (< i 256)
+              (aset syms i (->SymbolState 0 0 0 0 0))
+              (let [f (aget freqs i)]
+                (when (> f 0)
+                  (init! (aget syms i) total f))
+                (recur (inc i) (+ total f)))))
+        freq-table-size (write-frequencies0 out freqs)
+        _ (.rewind in)
+        compressed-data-size (encode-payload0 in syms out)]
+    (+ freq-table-size compressed-data-size)))
+
+(defn- encode1 ^long [^Buffer in ^ByteBuffer out]
+  (let [freqs (calculate-frequencies1 in)
+        ^objects syms (make-array SymbolState 256 256)
+        _ (dotimes [i 256]
+            (let [^ints fs (aget freqs i)]
+              (loop [j 0, total 0]
+                (when (< j 256)
+                  (aset ^objects (aget syms i) j (->SymbolState 0 0 0 0 0))
+                  (let [f (aget fs j)]
+                    (when (> f 0)
+                      (init! (aget2 syms i j) total f))
+                    (recur (inc j) (+ total f)))))))
+        freq-table-size (write-frequencies1 out freqs)
+        _ (.rewind in)
+        compressed-data-size (encode-payload1 in syms out)]
+    (+ freq-table-size compressed-data-size)))
+
+(defn encode
+  "Reads a byte sequence from the given ByteBuffer and encodes it by the rANS4x8 codec.
+  Returns the encoded result as a byte array."
+  ^bytes [^long order ^ByteBuffer in]
+  (let [raw-size (.remaining in)
+        ^ByteBuffer out (doto ^Buffer (allocate-output-buffer raw-size)
+                          (.mark)
+                          (.position PREFIX_BYTE_LEN))
+        ;; According to the specification, Order-1 encoding cannot be applicable
+        ;; to the input that is smaller than 4 bytes. So, in that case, the encoder
+        ;; automatically falls back to Order-0 encoding.
+        order' (if (< raw-size 4) 0 order)
+        compressed-size (case order'
+                          0 (encode0 in out)
+                          1 (encode1 in out))]
+    (.reset ^Buffer out)
+    (.put out (unchecked-byte order'))
+    (.putInt out compressed-size)
+    (.putInt out raw-size)
+    (Arrays/copyOfRange (.array out) 0 (+ PREFIX_BYTE_LEN compressed-size))))
diff --git a/src/cljam/io/cram/data_series.clj b/src/cljam/io/cram/data_series.clj
index 8ed6515..4ffc949 100644
--- a/src/cljam/io/cram/data_series.clj
+++ b/src/cljam/io/cram/data_series.clj
@@ -151,17 +151,17 @@
    ;; except for the value for :compressor, which is used when compressing embedded
    ;; reference blocks
    :embedded-ref {:content-id 1, :codec :external, :compressor :gzip}
-   :BF {:content-id  2, :codec :external, :compressor :gzip}
-   :CF {:content-id  3, :codec :external, :compressor :gzip}
-   :RI {:content-id  4, :codec :external, :compressor :gzip}
-   :RL {:content-id  5, :codec :external, :compressor :gzip}
-   :AP {:content-id  6, :codec :external, :compressor :gzip}
-   :RG {:content-id  7, :codec :external, :compressor :gzip}
+   :BF {:content-id  2, :codec :external, :compressor :r4x8-o1}
+   :CF {:content-id  3, :codec :external, :compressor :r4x8-o1}
+   :RI {:content-id  4, :codec :external, :compressor :r4x8-o0}
+   :RL {:content-id  5, :codec :external, :compressor :r4x8-o1}
+   :AP {:content-id  6, :codec :external, :compressor :r4x8-o0}
+   :RG {:content-id  7, :codec :external, :compressor :r4x8-o1}
    :RN {:content-id  8, :codec :byte-array-stop, :stop-byte (int \tab), :compressor :gzip}
    :MF {:content-id  9, :codec :external, :compressor :gzip}
-   :NS {:content-id 10, :codec :external, :compressor :gzip}
+   :NS {:content-id 10, :codec :external, :compressor :r4x8-o1}
    :NP {:content-id 11, :codec :external, :compressor :gzip}
-   :TS {:content-id 12, :codec :external, :compressor :gzip}
+   :TS {:content-id 12, :codec :external, :compressor :r4x8-o1}
    :NF {:content-id 13, :codec :external, :compressor :gzip}
    :TL {:content-id 14, :codec :external, :compressor :gzip}
    :FN {:content-id 15, :codec :external, :compressor :gzip}
@@ -185,8 +185,8 @@
         :len-encoding {:codec :external, :content-id 29, :compressor :gzip}
         :val-encoding {:codec :external, :content-id 30, :compressor :gzip}}
    :MQ {:content-id 31, :codec :external, :compressor :gzip}
-   :BA {:content-id 32, :codec :external, :compressor :gzip}
-   :QS {:content-id 33, :codec :external, :compressor :gzip}})
+   :BA {:content-id 32, :codec :external, :compressor :r4x8-o1}
+   :QS {:content-id 33, :codec :external, :compressor :r4x8-o1}})
 
 (defn- build-codec-encoder
   [{:keys [codec content-id compressor] :as params} data-type content-id->state]
diff --git a/src/cljam/io/cram/encode/compressor.clj b/src/cljam/io/cram/encode/compressor.clj
index d36a1fd..0e8a08d 100644
--- a/src/cljam/io/cram/encode/compressor.clj
+++ b/src/cljam/io/cram/encode/compressor.clj
@@ -1,4 +1,6 @@
 (ns cljam.io.cram.encode.compressor
+  (:require [cljam.io.cram.codecs.rans4x8 :as rans]
+            [cljam.io.util.byte-buffer :as bb])
   (:import [java.io ByteArrayOutputStream OutputStream]
            [org.apache.commons.compress.compressors.bzip2 BZip2CompressorOutputStream]
            [org.apache.commons.compress.compressors.gzip GzipCompressorOutputStream]
@@ -39,6 +41,16 @@
     (.finish out)
     {:compressor :lzma, :data (.toByteArray baos)}))
 
+(defn- compress-with-rans4x8 [order uncompressed]
+  (rans/encode order (bb/make-lsb-byte-buffer uncompressed)))
+
+(deftype RANS4x8Compressor [^long order ^ByteArrayOutputStream out]
+  ICompressor
+  (compressor-output-stream [_] out)
+  (->compressed-result [_]
+    {:compressor (if (zero? order) :r4x8-o0 :r4x8-o1),
+     :data (compress-with-rans4x8 order (.toByteArray out))}))
+
 (defn- compress-with [f ^bytes uncompressed]
   (let [out (ByteArrayOutputStream.)]
     (with-open [^OutputStream os (f out)]
@@ -62,6 +74,8 @@
                                                           uncompressed)
                                      :lzma (compress-with #(XZCompressorOutputStream. %)
                                                           uncompressed)
+                                     :r4x8-o0 (compress-with-rans4x8 0 uncompressed)
+                                     :r4x8-o1 (compress-with-rans4x8 1 uncompressed)
                                      (throw
                                       (ex-info (str "compression method " method
                                                     " not supported")
@@ -79,7 +93,9 @@
                 :gzip (->GzipCompressor (GzipCompressorOutputStream. baos) baos)
                 :bzip (->BZip2Compressor (BZip2CompressorOutputStream. baos) baos)
                 :lzma (->LZMACompressor (XZCompressorOutputStream. baos) baos)
-                :best (->SelectiveCompressor baos #{:raw :gzip :bzip :lzma})
+                :r4x8-o0 (->RANS4x8Compressor 0 baos)
+                :r4x8-o1 (->RANS4x8Compressor 1 baos)
+                :best (->SelectiveCompressor baos #{:raw :gzip :bzip :lzma :r4x8-o0 :r4x8-o1})
                 (if (set? method-or-methods)
                   (->SelectiveCompressor baos method-or-methods)
                   (throw
diff --git a/src/cljam/io/cram/encode/structure.clj b/src/cljam/io/cram/encode/structure.clj
index c00a05c..a52332f 100644
--- a/src/cljam/io/cram/encode/structure.clj
+++ b/src/cljam/io/cram/encode/structure.clj
@@ -95,7 +95,7 @@
 (defn encode-block
   "Encodes a block to the given OutputStream."
   [^OutputStream out method content-type content-id raw-size ^bytes block-data]
-  (let [method' (case method :raw 0 :gzip 1 :bzip 2 :lzma 3)]
+  (let [method' (case method :raw 0 :gzip 1 :bzip 2 :lzma 3 (:r4x8-o0 :r4x8-o1) 4)]
     (with-crc-suffixed out
       (fn [out']
         (lsb/write-ubyte out' method')


Test patch:
diff --git a/test/cljam/io/cram/codecs/rans4x8_test.clj b/test/cljam/io/cram/codecs/rans4x8_test.clj
index 9d83c16..466d47e 100644
--- a/test/cljam/io/cram/codecs/rans4x8_test.clj
+++ b/test/cljam/io/cram/codecs/rans4x8_test.clj
@@ -7,6 +7,9 @@
            [java.nio ByteBuffer]
            [java.util Arrays]))
 
+(defn- arr->freq-map [arr]
+  (into {} (keep-indexed (fn [i x] (when (pos? x) [i x]))) arr))
+
 (deftest read-frequencies0-test
   ;; cited from https://github.com/samtools/hts-specs/blob/346a94a9980b0105e926a019b4e62fa1b9e30910/CRAMcodecs.tex#L389-L394
   (let [arr (byte-array [0x61      0x87 0x47
@@ -21,10 +24,23 @@
             0x63 0x174
             0x64 0x174
             0x72 0x2e8}
-           (into {} (keep-indexed (fn [i x] (when (pos? x) [i x])))
-                 (#'rans/read-frequencies0 bb))))
+           (arr->freq-map (#'rans/read-frequencies0 bb))))
     (is (zero? (.remaining bb)))))
 
+(defn- arr->freq-map2 [arr]
+  (into {}
+        (keep-indexed
+         (fn [i arr]
+           (when-let [m (->> arr
+                             (into {}
+                                   (keep-indexed
+                                    (fn [j x]
+                                      (when (pos? x)
+                                        [j x]))))
+                             not-empty)]
+             [i m])))
+        arr))
+
 (deftest read-frequencies1-test
   ;; cited from https://github.com/samtools/hts-specs/blob/346a94a9980b0105e926a019b4e62fa1b9e30910/CRAMcodecs.tex#L489-L516
   (let [arr (byte-array [0x00
@@ -63,18 +79,7 @@
             0x63 {0x61 0xfff}
             0x64 {0x61 0xfff}
             0x72 {0x61 0xfff}}
-           (into {}
-                 (keep-indexed
-                  (fn [i arr]
-                    (when-let [m (->> arr
-                                      (into {}
-                                            (keep-indexed
-                                             (fn [j x]
-                                               (when (pos? x)
-                                                 [j x]))))
-                                      not-empty)]
-                      [i m])))
-                 (#'rans/read-frequencies1 bb))))
+           (arr->freq-map2 (#'rans/read-frequencies1 bb))))
     (is (zero? (.remaining bb)))))
 
 (deftest cumulative-frequencies-test
@@ -122,3 +127,104 @@
           compressed (read-as-buffer "test-resources/cram/codecs/rans4x8/compressed_order1.dat.rans4x8")]
       (is (Arrays/equals ^bytes (.array uncompressed) ^bytes (rans/decode compressed)))
       (is (zero? (.remaining compressed))))))
+
+(deftest calculate-frequencies0-test
+  (let [arr (.getBytes "abracadabra")
+        bb (bb/make-lsb-byte-buffer arr)]
+    (is (= {0x61 0x747
+            0x62 0x2e8
+            0x63 0x174
+            0x64 0x174
+            0x72 0x2e8}
+           (arr->freq-map (#'rans/calculate-frequencies0 bb)))))
+  (let [arr (->> (concat (repeat 6144 1)
+                         (repeat 2048 2)
+                         (repeat    1 3))
+                 shuffle
+                 byte-array)
+        bb (bb/make-lsb-byte-buffer arr)]
+    (is (= {1 3070
+            2 1024
+            3 1}
+           (arr->freq-map (#'rans/calculate-frequencies0 bb))))))
+
+(deftest calculate-frequencies1-test
+  (let [arr (.getBytes "abracadabraabracadabraabracadabraabracadabr")
+        bb (bb/make-lsb-byte-buffer arr)]
+    (is (= {0x00 {0x61 0x7ff
+                  0x62 0x400
+                  0x72 0x400}
+            0x61 {0x61 0x286
+                  0x62 0x6bd
+                  0x63 0x35e
+                  0x64 0x35e}
+            0x62 {0x72 0xfff}
+            0x63 {0x61 0xfff}
+            0x64 {0x61 0xfff}
+            0x72 {0x61 0xfff}}
+           (arr->freq-map2 (#'rans/calculate-frequencies1 bb)))))
+  (let [arr (->> (concat (repeat 6144 1)
+                         (repeat 2048 2)
+                         (repeat    1 3))
+                 shuffle
+                 (interleave (repeat 100))
+                 byte-array)
+        bb (bb/make-lsb-byte-buffer arr)]
+    (is (= {100 {1 3070
+                 2 1024
+                 3 1}
+            1 {100 4095}
+            2 {100 4095}
+            3 {100 4095}}
+           (-> (arr->freq-map2 (#'rans/calculate-frequencies1 bb))
+               (dissoc 0))))))
+
+(deftest write-frequencies0-test
+  (let [freqs {0x61 1863
+               0x62 744
+               0x63 372
+               0x64 372
+               0x72 744}
+        arr (int-array 256)
+        _ (doseq [[s f] freqs]
+            (aset arr s f))
+        bb (#'rans/allocate-output-buffer 11)]
+    (#'rans/write-frequencies0 bb arr)
+    (.flip bb)
+    (is (= freqs (arr->freq-map (#'rans/read-frequencies0 bb))))))
+
+(deftest write-frequencies1-test
+  (let [freqs {0x00 {0x61 0x7ff
+                     0x62 0x400
+                     0x72 0x400}
+               0x61 {0x61 0x286
+                     0x62 0x6bd
+                     0x63 0x35e
+                     0x64 0x35e}
+               0x62 {0x72 0xfff}
+               0x63 {0x61 0xfff}
+               0x64 {0x61 0xfff}
+               0x72 {0x61 0xfff}}
+        ^"[[I" arr (make-array Integer/TYPE 256 256)
+        _ (doseq [[s1 fs] freqs
+                  [s2 f] fs]
+            (aset ^ints (aget arr s1) s2 f))
+        bb (#'rans/allocate-output-buffer 43)]
+    (#'rans/write-frequencies1 bb arr)
+    (.flip bb)
+    (is (= freqs (arr->freq-map2 (#'rans/read-frequencies1 bb))))))
+
+(deftest encode-test
+  (testing "Order-0"
+    (let [uncompressed (read-as-buffer "test-resources/cram/codecs/rans4x8/uncompressed_order0.dat")
+          compressed (bb/make-lsb-byte-buffer (rans/encode 0 uncompressed))]
+      (is (Arrays/equals (.array uncompressed) (rans/decode compressed)))))
+  (testing "Order-1"
+    (let [uncompressed (read-as-buffer "test-resources/cram/codecs/rans4x8/uncompressed_order1.dat")
+          compressed (bb/make-lsb-byte-buffer (rans/encode 1 uncompressed))]
+      (is (Arrays/equals (.array uncompressed) (rans/decode compressed)))))
+  (testing "Specifying Order-1, but ended up falling back to Order-0"
+    (let [uncompressed (bb/make-lsb-byte-buffer (byte-array (.getBytes "ab")))
+          compressed (bb/make-lsb-byte-buffer (rans/encode 1 uncompressed))]
+      (is (.get compressed 0))
+      (is (Arrays/equals (.array uncompressed) (rans/decode compressed))))))