datasketches 0.5.0

// Licensed to the Apache Software Foundation (ASF) under one
// or more contributor license agreements.  See the NOTICE file
// distributed with this work for additional information
// regarding copyright ownership.  The ASF licenses this file
// to you under the Apache License, Version 2.0 (the
// "License"); you may not use this file except in compliance
// with the License.  You may obtain a copy of the License at
//
//   http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing,
// software distributed under the License is distributed on an
// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
// KIND, either express or implied.  See the License for the
// specific language governing permissions and limitations
// under the License.

use std::hash::Hash;

use crate::codec::SketchBytes;
use crate::codec::SketchSlice;
use crate::codec::assert::ensure_preamble_longs_in;
use crate::codec::assert::ensure_serial_version_is;
use crate::codec::assert::insufficient_data;
use crate::codec::family::Family;
use crate::common::NumStdDev;
use crate::common::inv_pow2::inv_pow2;
use crate::cpc::DEFAULT_LG_K;
use crate::cpc::Flavor;
use crate::cpc::MAX_LG_K;
use crate::cpc::MIN_LG_K;
use crate::cpc::compression::decode_pairs;
use crate::cpc::compression::decode_window;
use crate::cpc::compression::determine_pseudo_phase;
use crate::cpc::compression::encode_pairs;
use crate::cpc::compression::encode_window;
use crate::cpc::compression_data::COLUMN_PERMUTATIONS_FOR_DECODING;
use crate::cpc::compression_data::COLUMN_PERMUTATIONS_FOR_ENCODING;
use crate::cpc::determine_correct_offset;
use crate::cpc::determine_flavor;
use crate::cpc::estimator::estimate;
use crate::cpc::estimator::lower_bound;
use crate::cpc::estimator::upper_bound;
use crate::cpc::kxp_byte_lookup::KXP_BYTE_TABLE;
use crate::cpc::pair_table::PairTable;
use crate::cpc::serialization::FLAG_COMPRESSED;
use crate::cpc::serialization::FLAG_HAS_HIP;
use crate::cpc::serialization::FLAG_HAS_TABLE;
use crate::cpc::serialization::FLAG_HAS_WINDOW;
use crate::cpc::serialization::SERIAL_VERSION;
use crate::cpc::serialization::make_preamble_ints;
use crate::error::Error;
use crate::error::ErrorKind;
use crate::hash::DEFAULT_UPDATE_SEED;
use crate::hash::MurmurHash3X64128;
use crate::hash::check_seed_hash;
use crate::hash::compute_seed_hash;

/// A Compressed Probabilistic Counting sketch.
///
/// See the [module level documentation](super) for more.
#[derive(Debug, Clone)]
pub struct CpcSketch {
    // immutable config variables
    lg_k: u8,
    seed: u64,
    seed_hash: u16,

    // sketch state
    /// Part of a speed optimization.
    pub(super) first_interesting_column: u8,
    /// The number of coupons collected so far.
    pub(super) num_coupons: u32,
    /// Sparse and surprising values.
    pub(super) surprising_value_table: Option<PairTable>,
    /// Derivable from num_coupons, but made explicit for speed.
    pub(super) window_offset: u8,
    /// Size K bytes in dense mode (flavor >= HYBRID).
    pub(super) sliding_window: Vec<u8>,

    // estimator state
    /// Whether the sketch is a result of merging.
    ///
    /// If `false`, the HIP (Historical Inverse Probability) estimator is used.
    /// If `true`, the ICON (Inter-Column Optimal) Estimator is fallback in use.
    pub(super) merge_flag: bool,
    // the following variables are only valid in HIP estimator
    /// A pre-calculated probability factor (`k * p`) used to compute the increment delta.
    kxp: f64,
    /// The accumulated cardinality estimate.
    hip_est_accum: f64,
}

impl Default for CpcSketch {
    fn default() -> Self {
        Self::new(DEFAULT_LG_K).unwrap()
    }
}

impl CpcSketch {
    /// Creates a new `CpcSketch` with the given `lg_k` and default seed.
    ///
    /// # Errors
    ///
    /// Returns an error if `lg_k` is not in the range `[4, 26]`.
    pub fn new(lg_k: u8) -> Result<Self, Error> {
        Self::with_seed(lg_k, DEFAULT_UPDATE_SEED)
    }

    /// Creates a new `CpcSketch` with the given `lg_k` and `seed`.
    ///
    /// # Errors
    ///
    /// Returns an error if `lg_k` is not in the range `[4, 26]`, or the computed seed hash is zero.
    pub fn with_seed(lg_k: u8, seed: u64) -> Result<Self, Error> {
        if !(MIN_LG_K..=MAX_LG_K).contains(&lg_k) {
            return Err(Error::invalid_argument(format!(
                "lg_k must be in [{MIN_LG_K}, {MAX_LG_K}], got {lg_k}"
            )));
        }

        Ok(Self {
            lg_k,
            seed,
            seed_hash: compute_seed_hash(seed, ErrorKind::InvalidArgument)?,
            first_interesting_column: 0,
            num_coupons: 0,
            surprising_value_table: None,
            window_offset: 0,
            sliding_window: vec![],
            merge_flag: false,
            kxp: (1 << lg_k) as f64,
            hip_est_accum: 0.0,
        })
    }

    /// Returns the configured `lg_k`.
    pub fn lg_k(&self) -> u8 {
        self.lg_k
    }

    /// Returns the best estimate of the cardinality of the sketch.
    pub fn estimate(&self) -> f64 {
        estimate(
            self.merge_flag,
            self.hip_est_accum,
            self.lg_k,
            self.num_coupons,
        )
    }

    /// Returns the best estimate of the lower bound of the confidence interval.
    pub fn lower_bound(&self, num_std_dev: NumStdDev) -> f64 {
        lower_bound(
            self.merge_flag,
            self.hip_est_accum,
            self.lg_k,
            self.num_coupons,
            num_std_dev,
        )
    }

    /// Returns the best estimate of the upper bound of the confidence interval.
    pub fn upper_bound(&self, num_std_dev: NumStdDev) -> f64 {
        upper_bound(
            self.merge_flag,
            self.hip_est_accum,
            self.lg_k,
            self.num_coupons,
            num_std_dev,
        )
    }

    /// Returns `true` if the sketch is empty.
    pub fn is_empty(&self) -> bool {
        self.num_coupons == 0
    }

    /// Updates the sketch with a hashable value.
    ///
    /// You may use [`hash::value`](crate::hash::value) wrappers when another DataSketches
    /// implementation requires a specific value hashing strategy.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use datasketches::cpc::CpcSketch;
    /// use datasketches::hash::value::canonical_float;
    ///
    /// let mut sketch = CpcSketch::with_seed(11, 123).unwrap();
    /// sketch.update(1);
    /// sketch.update(2);
    /// sketch.update(3);
    ///
    /// let mut sketch = CpcSketch::with_seed(11, 123).unwrap();
    /// sketch.update(canonical_float::from_f64(1.5));
    /// sketch.update(canonical_float::from_f64(2.5));
    /// sketch.update(canonical_float::from_f64(3.5));
    /// ```
    pub fn update<T: Hash>(&mut self, value: T) {
        let mut hasher = MurmurHash3X64128::with_seed(self.seed);
        value.hash(&mut hasher);
        let (h1, h2) = hasher.finish128();

        let k = 1 << self.lg_k;
        let col = h2.leading_zeros(); // 0 <= col <= 64
        let col = if col > 63 { 63 } else { col as u8 }; // clip so that 0 <= col <= 63
        let row = (h1 & (k - 1)) as u32;
        let mut row_col = (row << 6) | (col as u32);
        // To avoid the hash table's "empty" value, we change the row of the following pair.
        // This case is extremely unlikely, but we might as well handle it.
        if row_col == u32::MAX {
            row_col ^= 1 << 6;
        }
        self.row_col_update(row_col);
    }

    pub(super) fn flavor(&self) -> Flavor {
        determine_flavor(self.lg_k, self.num_coupons)
    }

    pub(super) fn row_col_update(&mut self, row_col: u32) {
        let col = (row_col & 63) as u8;
        if col < self.first_interesting_column {
            // important speed optimization
            return;
        }

        if self.num_coupons == 0 {
            // promote EMPTY to SPARSE
            self.surprising_value_table = Some(PairTable::new(2, 6 + self.lg_k));
        }

        if self.sliding_window.is_empty() {
            self.update_sparse(row_col);
        } else {
            self.update_windowed(row_col);
        }
    }

    pub(super) fn seed(&self) -> u64 {
        self.seed
    }

    pub(super) fn surprising_value_table(&self) -> &PairTable {
        self.surprising_value_table
            .as_ref()
            .expect("surprising value table must be initialized")
    }

    fn surprising_value_table_mut(&mut self) -> &mut PairTable {
        self.surprising_value_table
            .as_mut()
            .expect("surprising value table must be initialized")
    }

    fn update_hip(&mut self, row_col: u32) {
        let k = 1 << self.lg_k;
        let col = (row_col & 63) as usize;
        let one_over_p = (k as f64) / self.kxp;
        self.hip_est_accum += one_over_p;
        self.kxp -= inv_pow2((col + 1) as u8) // notice the "+1"
    }

    fn update_sparse(&mut self, row_col: u32) {
        let k = 1 << self.lg_k;
        let c32pre = (self.num_coupons as u64) << 5;
        debug_assert!(c32pre < 3 * k); // C < 3K/32, in other words, flavor == SPARSE
        let is_novel = self.surprising_value_table_mut().maybe_insert(row_col);
        if is_novel {
            self.num_coupons += 1;
            self.update_hip(row_col);
            let c32post = (self.num_coupons as u64) << 5;
            if c32post >= 3 * k {
                self.promote_sparse_to_windowed();
            }
        }
    }

    fn promote_sparse_to_windowed(&mut self) {
        debug_assert_eq!(self.window_offset, 0);

        let k = 1 << self.lg_k;
        let c32 = (self.num_coupons as u64) << 5;
        debug_assert!((c32 == (3 * k)) || ((self.lg_k == 4) && (c32 > (3 * k))));

        self.sliding_window.resize(k as usize, 0);

        let old_table = self
            .surprising_value_table
            .replace(PairTable::new(2, 6 + self.lg_k))
            .expect("surprising value table must be initialized");
        let old_slots = old_table.slots();
        for &row_col in old_slots {
            if row_col != u32::MAX {
                let col = (row_col & 63) as u8;
                if col < 8 {
                    let row = (row_col >> 6) as usize;
                    self.sliding_window[row] |= 1 << col;
                } else {
                    // cannot use must_insert(), because it doesn't provide for growth
                    let is_novel = self.surprising_value_table_mut().maybe_insert(row_col);
                    debug_assert!(is_novel);
                }
            }
        }
    }

    fn update_windowed(&mut self, row_col: u32) {
        debug_assert!(self.window_offset <= 56);
        let k = 1 << self.lg_k;
        let c32pre = (self.num_coupons as u64) << 5;
        debug_assert!(c32pre >= 3 * k); // C >= 3K/32, in other words flavor >= HYBRID
        let c8pre = (self.num_coupons as u64) << 3;
        let w8pre = (self.window_offset as u64) << 3;
        debug_assert!(c8pre < (27 + w8pre) * k); // C < (K * 27/8) + (K * windowOffset)

        let mut is_novel = false; // novel if new coupon;
        let col = (row_col & 63) as u8;
        if col < self.window_offset {
            // track the surprising 0's "before" the window
            is_novel = self.surprising_value_table_mut().maybe_delete(row_col); // inverted logic
        } else if col < self.window_offset + 8 {
            // track the 8 bits inside the window
            let row = (row_col >> 6) as usize;
            let old_bits = self.sliding_window[row];
            let new_bits = old_bits | (1 << (col - self.window_offset));
            if old_bits != new_bits {
                self.sliding_window[row] = new_bits;
                is_novel = true;
            }
        } else {
            // track the surprising 1's "after" the window
            is_novel = self.surprising_value_table_mut().maybe_insert(row_col); // normal logic
        }

        if is_novel {
            self.num_coupons += 1;
            self.update_hip(row_col);
            let c8post = (self.num_coupons as u64) << 3;
            if c8post >= (27 + w8pre) * k {
                self.move_window();
                debug_assert!((1..=56).contains(&self.window_offset));
                let w8post = (self.window_offset as u64) << 3;
                debug_assert!(c8post < ((27 + w8post) * k)); // C < (K * 27/8) + (K * windowOffset)
            }
        }
    }

    fn move_window(&mut self) {
        let new_offset = self.window_offset + 1;
        debug_assert!(new_offset <= 56);
        debug_assert_eq!(
            new_offset,
            determine_correct_offset(self.lg_k, self.num_coupons)
        );

        let k = 1 << self.lg_k;

        // Construct the full-sized bit matrix that corresponds to the sketch
        let bit_matrix = self.build_bit_matrix();

        // refresh the KXP register on every 8th window shift.
        if (new_offset & 0x7) == 0 {
            self.refresh_kxp(&bit_matrix);
        }

        self.surprising_value_table_mut().clear(); // the new number of surprises will be about the same

        let mask_for_clearing_window = (0xFF << new_offset) ^ u64::MAX;
        let mask_for_flipping_early_zone = (1u64 << new_offset) - 1;

        let mut all_surprises_ored = 0u64;
        for i in 0..k {
            let mut pattern = bit_matrix[i];
            self.sliding_window[i] = ((pattern >> new_offset) & 0xff) as u8;
            pattern &= mask_for_clearing_window;
            // The following line converts surprising 0's to 1's in the "early zone",
            // (and vice versa, which is essential for this procedure's O(k) time cost).
            pattern ^= mask_for_flipping_early_zone;
            all_surprises_ored |= pattern; // a cheap way to recalculate first_interesting_column
            while pattern != 0 {
                let col = pattern.trailing_zeros();
                pattern ^= 1 << col; // erase the 1
                let row_col = ((i as u32) << 6) | col;
                let is_novel = self.surprising_value_table_mut().maybe_insert(row_col);
                debug_assert!(is_novel);
            }
        }

        self.window_offset = new_offset;
        self.first_interesting_column = all_surprises_ored.trailing_zeros() as u8;
        if self.first_interesting_column > new_offset {
            self.first_interesting_column = new_offset; // corner case
        }
    }

    /// The KXP register is a double with roughly 50 bits of precision, but it might need roughly
    /// 90 bits to track the value with perfect accuracy.
    ///
    /// Therefore, we recalculate KXP occasionally from the sketch's full bit_matrix so that it
    /// will reflect changes that were previously outside the mantissa.
    fn refresh_kxp(&mut self, bit_matrix: &[u64]) {
        // for improved numerical accuracy, we separately sum the bytes of the u64's
        let mut byte_sums = [0.0; 8];
        for &bit in bit_matrix {
            let mut word = bit;
            for sum in byte_sums.iter_mut() {
                let byte = (word & 0xFF) as usize;
                *sum += KXP_BYTE_TABLE[byte];
                word >>= 8;
            }
        }

        let mut total = 0.0;
        for i in (0..8).rev() {
            // the reverse order is important
            let factor = inv_pow2((i * 8) as u8); // pow (256.0, (-1.0 * ((double) j)))
            total += factor * byte_sums[i];
        }

        self.kxp = total;
    }

    pub(super) fn build_bit_matrix(&self) -> Vec<u64> {
        let k = 1 << self.lg_k;
        let offset = self.window_offset;
        debug_assert!(offset <= 56);

        // Fill the matrix with default rows in which the "early zone" is filled with ones.
        // This is essential for the routine's O(k) time cost (as opposed to O(C)).
        let default_row = (1u64 << offset) - 1;

        let mut matrix = vec![default_row; k];
        if self.num_coupons == 0 {
            return matrix;
        }

        if !self.sliding_window.is_empty() {
            // In other words, we are in window mode, not sparse mode
            for i in 0..k {
                // set the window bits, trusting the sketch's current offset
                matrix[i] |= (self.sliding_window[i] as u64) << offset;
            }
        }

        for &row_col in self.surprising_value_table().slots() {
            if row_col != u32::MAX {
                let col = (row_col & 63) as u8;
                let row = (row_col >> 6) as usize;
                // Flip the specified matrix bit from its default value.
                // In the "early" zone the bit changes from 1 to 0.
                // In the "late" zone the bit changes from 0 to 1.
                matrix[row] ^= 1 << col;
            }
        }

        matrix
    }

    /// Returns the estimated size of the sketch in bytes.
    pub fn estimated_size(&self) -> usize {
        let heap_size = self.sliding_window.capacity()
            + self
                .surprising_value_table
                .as_ref()
                .map(|t| t.estimated_size())
                .unwrap_or(0);

        size_of::<Self>() + heap_size
    }
}

impl CpcSketch {
    /// Serializes this `CpcSketch` to bytes.
    pub fn serialize(&self) -> Vec<u8> {
        let flavor = self.flavor();
        let has_hip = !self.merge_flag;
        let has_window = matches!(flavor, Flavor::Pinned | Flavor::Sliding);
        let mut pairs = match flavor {
            Flavor::Empty => vec![],
            Flavor::Sparse => {
                debug_assert!(self.sliding_window.is_empty());
                self.surprising_value_table().unwrapping_get_items()
            }
            Flavor::Hybrid => {
                debug_assert!(!self.sliding_window.is_empty());
                debug_assert_eq!(self.window_offset, 0);

                let mut table_pairs = self.surprising_value_table().unwrapping_get_items();
                table_pairs.sort_unstable();
                let num_table_pairs = table_pairs.len();
                let mut all_pairs = vec![0; self.num_coupons as usize];

                let mut index = num_table_pairs;
                for (row_index, &window_byte) in self.sliding_window.iter().enumerate() {
                    let mut window_byte = window_byte;
                    while window_byte != 0 {
                        let col_index = window_byte.trailing_zeros();
                        window_byte ^= 1 << col_index;
                        all_pairs[index] = ((row_index << 6) as u32) | col_index;
                        index += 1;
                    }
                }
                assert_eq!(index, all_pairs.len());

                let mut table_index = 0;
                let mut window_index = num_table_pairs;
                for final_index in 0..all_pairs.len() {
                    if table_index < num_table_pairs
                        && (window_index >= all_pairs.len()
                            || table_pairs[table_index] <= all_pairs[window_index])
                    {
                        all_pairs[final_index] = table_pairs[table_index];
                        table_index += 1;
                    } else {
                        all_pairs[final_index] = all_pairs[window_index];
                        window_index += 1;
                    }
                }
                all_pairs
            }
            Flavor::Pinned => {
                let mut pairs = self.surprising_value_table().unwrapping_get_items();
                for pair in &mut pairs {
                    assert!(*pair & 63 >= 8, "pair column index is less than 8: {pair}");
                    *pair -= 8;
                }
                pairs
            }
            Flavor::Sliding => {
                let mut pairs = self.surprising_value_table().unwrapping_get_items();
                let pseudo_phase = determine_pseudo_phase(self.lg_k, self.num_coupons);
                let permutation = &COLUMN_PERMUTATIONS_FOR_ENCODING[pseudo_phase as usize];
                debug_assert!(self.window_offset <= 56);
                for pair in &mut pairs {
                    let row = *pair >> 6;
                    let col = ((*pair & 63) as u8 + 56 - self.window_offset) & 63;
                    debug_assert!(col < 56);
                    *pair = (row << 6) | u32::from(permutation[col as usize]);
                }
                pairs
            }
        };
        pairs.sort_unstable();

        let table_num_entries = pairs.len() as u32;
        let mut payload = SketchBytes::with_capacity(if self.is_empty() { 0 } else { 256 });
        let window_words = has_window.then(|| {
            encode_window(
                &self.sliding_window,
                self.lg_k,
                self.num_coupons,
                &mut payload,
            )
        });
        let table_words =
            (!pairs.is_empty()).then(|| encode_pairs(&pairs, self.lg_k, &mut payload));
        let payload = payload.into_bytes();

        let has_table = table_words.is_some();
        let preamble_ints = make_preamble_ints(self.num_coupons, has_hip, has_table, has_window);
        let mut bytes = SketchBytes::with_capacity(40 + payload.len());
        bytes.write_u8(preamble_ints);
        bytes.write_u8(SERIAL_VERSION);
        bytes.write_u8(Family::CPC.id);
        bytes.write_u8(self.lg_k);
        bytes.write_u8(self.first_interesting_column);
        let flags = (1 << FLAG_COMPRESSED)
            | (if has_hip { 1 } else { 0 } << FLAG_HAS_HIP)
            | (if has_table { 1 } else { 0 } << FLAG_HAS_TABLE)
            | (if has_window { 1 } else { 0 } << FLAG_HAS_WINDOW);
        bytes.write_u8(flags);
        debug_assert_eq!(
            self.seed_hash,
            compute_seed_hash(self.seed, ErrorKind::InvalidArgument).unwrap()
        );
        bytes.write_u16_le(self.seed_hash);
        if !self.is_empty() {
            bytes.write_u32_le(self.num_coupons);
            if has_table && has_window {
                // if there is no window it is the same as number of coupons
                bytes.write_u32_le(table_num_entries);
                // HIP values can be in two different places in the sequence of fields
                // this is the first HIP decision point
                if has_hip {
                    self.write_hip(&mut bytes);
                }
            }
            if let Some(table_words) = table_words {
                debug_assert!(table_words <= u32::MAX as usize);
                bytes.write_u32_le(table_words as u32);
            }
            if let Some(window_words) = window_words {
                debug_assert!(window_words <= u32::MAX as usize);
                bytes.write_u32_le(window_words as u32);
            }
            // this is the second HIP decision point
            if has_hip && !(has_table && has_window) {
                self.write_hip(&mut bytes);
            }
            bytes.write(&payload);
        }
        bytes.into_bytes()
    }

    /// Deserializes a `CpcSketch` from bytes.
    ///
    /// # Errors
    ///
    /// Returns `InvalidData` if the image is malformed or its seed hash does not match the default
    /// seed.
    pub fn deserialize(bytes: &[u8]) -> Result<Self, Error> {
        Self::deserialize_with_seed(bytes, DEFAULT_UPDATE_SEED)
    }

    /// Deserializes a `CpcSketch` from bytes with the provided seed.
    ///
    /// # Errors
    ///
    /// Returns `InvalidData` if the image is malformed, its seed hash does not match `seed`, or
    /// `seed` itself computes to the reserved zero seed hash.
    pub fn deserialize_with_seed(bytes: &[u8], seed: u64) -> Result<Self, Error> {
        let mut cursor = SketchSlice::new(bytes);
        let preamble_ints = cursor
            .read_u8()
            .map_err(insufficient_data("preamble_ints"))?;
        let serial_version = cursor
            .read_u8()
            .map_err(insufficient_data("serial_version"))?;
        let family_id = cursor.read_u8().map_err(insufficient_data("family_id"))?;
        Family::CPC.validate_id(family_id)?;
        ensure_serial_version_is(SERIAL_VERSION, serial_version)?;

        let lg_k = cursor.read_u8().map_err(insufficient_data("lg_k"))?;
        let first_interesting_column = cursor
            .read_u8()
            .map_err(insufficient_data("first_interesting_column"))?;

        let flags = cursor.read_u8().map_err(insufficient_data("flags"))?;
        let seed_hash = cursor
            .read_u16_le()
            .map_err(insufficient_data("seed_hash"))?;
        let is_compressed = flags & (1 << FLAG_COMPRESSED) != 0;
        if !is_compressed {
            return Err(Error::deserial("only compressed sketches are supported"));
        }
        let has_hip = flags & (1 << FLAG_HAS_HIP) != 0;
        let has_table = flags & (1 << FLAG_HAS_TABLE) != 0;
        let has_window = flags & (1 << FLAG_HAS_WINDOW) != 0;

        let mut num_coupons = 0;
        let mut table_num_entries = 0;
        let mut table_data_words = 0;
        let mut window_data_words = 0;
        let mut kxp = 0.0;
        let mut hip_est_accum = 0.0;

        if has_table || has_window {
            num_coupons = cursor
                .read_u32_le()
                .map_err(insufficient_data("num_coupons"))?;
            if has_table && has_window {
                table_num_entries = cursor
                    .read_u32_le()
                    .map_err(insufficient_data("table_num_entries"))?;
                if has_hip {
                    kxp = cursor.read_f64_le().map_err(insufficient_data("kxp"))?;
                    hip_est_accum = cursor
                        .read_f64_le()
                        .map_err(insufficient_data("hip_est_accum"))?;
                }
            }
            if has_table {
                table_data_words = cursor
                    .read_u32_le()
                    .map_err(insufficient_data("table_data_words"))?
                    as usize;
            }
            if has_window {
                window_data_words = cursor
                    .read_u32_le()
                    .map_err(insufficient_data("window_data_words"))?
                    as usize;
            }
            if has_hip && !(has_table && has_window) {
                kxp = cursor.read_f64_le().map_err(insufficient_data("kxp"))?;
                hip_est_accum = cursor
                    .read_f64_le()
                    .map_err(insufficient_data("hip_est_accum"))?;
            }
            if !has_window {
                table_num_entries = num_coupons;
            }
        }

        let expected_preamble_ints =
            make_preamble_ints(num_coupons, has_hip, has_table, has_window);
        ensure_preamble_longs_in(&[expected_preamble_ints], preamble_ints)?;
        check_seed_hash(
            compute_seed_hash(seed, ErrorKind::InvalidData)?,
            seed_hash,
            "deserialized CpcSketch",
            ErrorKind::InvalidData,
        )?;
        if !(MIN_LG_K..=MAX_LG_K).contains(&lg_k) {
            return Err(Error::deserial(format!("lg_k out of range; got {}", lg_k)));
        }
        if first_interesting_column > 63 {
            return Err(Error::deserial(format!(
                "first_interesting_column out of range; got {}",
                first_interesting_column
            )));
        }

        // The coupon space of a sketch has `k * 64` cells (`k` rows of 64 columns each), so a
        // valid sketch can never report more coupons than that. Rejecting larger values keeps the
        // flavor arithmetic below from overflowing on corrupt input.
        if (num_coupons as u64) > 64 * (1u64 << lg_k) {
            return Err(Error::deserial(format!(
                "num_coupons ({}) exceeds coupon space for lg_k = {}",
                num_coupons, lg_k
            )));
        }

        // A valid sketch stores a sliding window exactly for the pinned and sliding flavors, and
        // stores its coupons in the surprising-value table for the sparse and hybrid flavors. The
        // flavor is fully determined by `lg_k` and `num_coupons`, so the flags must agree with it.
        let flavor = determine_flavor(lg_k, num_coupons);
        let window_expected = matches!(flavor, Flavor::Pinned | Flavor::Sliding);
        if has_window != window_expected {
            return Err(Error::deserial(format!(
                "sliding-window flag ({}) is inconsistent with the {:?} flavor",
                has_window, flavor
            )));
        }
        if matches!(flavor, Flavor::Sparse | Flavor::Hybrid) && !has_table {
            return Err(Error::deserial(format!(
                "table flag is unset but required for the {:?} flavor",
                flavor
            )));
        }

        // The number of stored table entries can never exceed the number of coupons.
        if table_num_entries > num_coupons {
            return Err(Error::deserial(format!(
                "table_num_entries ({}) exceeds num_coupons ({})",
                table_num_entries, num_coupons
            )));
        }
        // A pair requires at least one bit to encode, so the declared number of table entries can
        // never exceed the number of bits available in the table data. This also bounds the size
        // of the allocation made while decoding, rejecting corrupt inputs that claim an enormous
        // entry count backed by only a few data words.
        if (table_num_entries as usize) > table_data_words.saturating_mul(32) {
            return Err(Error::deserial(format!(
                "table_num_entries ({}) exceeds capacity of table data ({} words)",
                table_num_entries, table_data_words
            )));
        }
        let k = 1 << lg_k;
        if has_window && window_data_words.saturating_mul(32) < k {
            return Err(Error::deserial(format!(
                "window data ({} words) is too short for lg_k = {lg_k}",
                window_data_words
            )));
        }

        let window_data_bytes = window_data_words.checked_mul(4).ok_or_else(|| {
            Error::deserial("CPC window data word count overflows payload length")
        })?;
        let table_data_bytes = table_data_words
            .checked_mul(4)
            .ok_or_else(|| Error::deserial("CPC table data word count overflows payload length"))?;
        let payload_bytes = window_data_bytes
            .checked_add(table_data_bytes)
            .ok_or_else(|| Error::deserial("CPC payload length overflows"))?;
        let payload = cursor
            .remaining()
            .get(..payload_bytes)
            .ok_or_else(|| Error::deserial("insufficient data for CPC compressed payload"))?;
        let (window_data, table_data) = payload.split_at(window_data_bytes);
        let (table, window) = match flavor {
            Flavor::Empty => (PairTable::new(2, lg_k + 6), vec![]),
            Flavor::Sparse => {
                debug_assert!(window_data.is_empty(), "window is not expected");
                let pairs = decode_pairs(table_data, table_num_entries, lg_k)?;
                (
                    PairTable::from_slots(lg_k, table_num_entries, pairs)?,
                    vec![],
                )
            }
            Flavor::Hybrid => {
                debug_assert!(window_data.is_empty(), "window is not expected");
                let mut pairs = decode_pairs(table_data, table_num_entries, lg_k)?;
                let mut window = vec![0u8; 1 << lg_k];
                let mut next_true_pair = 0;
                for index in 0..table_num_entries as usize {
                    let row_col = pairs[index];
                    let col = row_col & 63;
                    if col < 8 {
                        window[(row_col >> 6) as usize] |= 1 << col;
                    } else {
                        pairs[next_true_pair as usize] = row_col;
                        next_true_pair += 1;
                    }
                }
                (PairTable::from_slots(lg_k, next_true_pair, pairs)?, window)
            }
            Flavor::Pinned => {
                let window = decode_window(window_data, lg_k, num_coupons)?;
                let table = if table_num_entries == 0 {
                    PairTable::new(2, lg_k + 6)
                } else {
                    let mut pairs = decode_pairs(table_data, table_num_entries, lg_k)?;
                    for pair in &mut pairs {
                        if (*pair & 63) >= 56 {
                            return Err(Error::deserial(format!(
                                "CPC pinned table pair column index is invalid: {pair}"
                            )));
                        }
                        *pair += 8;
                    }
                    PairTable::from_slots(lg_k, table_num_entries, pairs)?
                };
                (table, window)
            }
            Flavor::Sliding => {
                let window = decode_window(window_data, lg_k, num_coupons)?;
                let table = if table_num_entries == 0 {
                    PairTable::new(2, lg_k + 6)
                } else {
                    let mut pairs = decode_pairs(table_data, table_num_entries, lg_k)?;
                    let pseudo_phase = determine_pseudo_phase(lg_k, num_coupons);
                    let permutation = &COLUMN_PERMUTATIONS_FOR_DECODING[pseudo_phase as usize];
                    let offset = determine_correct_offset(lg_k, num_coupons);
                    if offset > 56 {
                        return Err(Error::deserial(format!(
                            "CPC sliding window offset is invalid: {offset}"
                        )));
                    }
                    for pair in &mut pairs {
                        let row = *pair >> 6;
                        let col = (*pair & 63) as usize;
                        if col >= permutation.len() {
                            return Err(Error::deserial(format!(
                                "CPC sliding table pair column index is invalid: {pair}"
                            )));
                        }
                        let col = (permutation[col] + offset + 8) & 63;
                        *pair = (row << 6) | u32::from(col);
                    }
                    PairTable::from_slots(lg_k, table_num_entries, pairs)?
                };
                (table, window)
            }
        };
        Ok(CpcSketch {
            lg_k,
            seed,
            seed_hash,
            first_interesting_column,
            num_coupons,
            surprising_value_table: Some(table),
            window_offset: determine_correct_offset(lg_k, num_coupons),
            sliding_window: window,
            merge_flag: !has_hip,
            kxp,
            hip_est_accum,
        })
    }

    fn write_hip(&self, bytes: &mut SketchBytes) {
        bytes.write_f64_le(self.kxp);
        bytes.write_f64_le(self.hip_est_accum);
    }
}

impl CpcSketch {
    /// Returns the estimated maximum compressed serialized size of a sketch.
    ///
    /// The actual size of a compressed CPC sketch has a small random variance, but the following
    /// empirically measured size should be large enough for at least 99.9 percent of sketches.
    ///
    /// For small values of `n` the size can be much smaller.
    ///
    /// # Errors
    ///
    /// Returns an error if `lg_k` is not in the range `[4, 26]`.
    pub fn max_serialized_bytes(lg_k: u8) -> Result<usize, Error> {
        if !(MIN_LG_K..=MAX_LG_K).contains(&lg_k) {
            return Err(Error::invalid_argument(format!(
                "lg_k must be in [{MIN_LG_K}, {MAX_LG_K}], got {lg_k}"
            )));
        }

        // These empirical values for the 99.9th percentile of size in bytes were measured using
        // 100,000 trials. The value for each trial is the maximum of 5*16=80 measurements
        // that were equally spaced over values of the quantity C/K between 3.0 and 8.0.
        // This table does not include the worst-case space for the preamble, which is added
        // by the function.
        const MAX_PREAMBLE_SIZE_BYTES: usize = 40;
        const EMPIRICAL_SIZE_MAX_LGK: u8 = 19;
        const EMPIRICAL_MAX_SIZE_FACTOR: f64 = 0.6; // 0.6 = 4.8 / 8.0
        static EMPIRICAL_MAX_SIZE_BYTES: [usize; 16] = [
            24,     // lg_k = 4
            36,     // lg_k = 5
            56,     // lg_k = 6
            100,    // lg_k = 7
            180,    // lg_k = 8
            344,    // lg_k = 9
            660,    // lg_k = 10
            1292,   // lg_k = 11
            2540,   // lg_k = 12
            5020,   // lg_k = 13
            9968,   // lg_k = 14
            19836,  // lg_k = 15
            39532,  // lg_k = 16
            78880,  // lg_k = 17
            157516, // lg_k = 18
            314656, // lg_k = 19
        ];

        let max_bytes = if lg_k <= EMPIRICAL_SIZE_MAX_LGK {
            EMPIRICAL_MAX_SIZE_BYTES[(lg_k - MIN_LG_K) as usize] + MAX_PREAMBLE_SIZE_BYTES
        } else {
            let k = 1 << lg_k;
            ((EMPIRICAL_MAX_SIZE_FACTOR * k as f64) as usize) + MAX_PREAMBLE_SIZE_BYTES
        };
        Ok(max_bytes)
    }
}