text: /**
   * Copyright (c) Facebook, Inc. and its affiliates.
   *
   * This source code is licensed under the MIT license found in the
   * LICENSE file in the root directory of this source tree.
   */
  
  // -*- c++ -*-
  
  #include <faiss/impl/AdditiveQuantizer.h>
  
  #include <cstddef>
  #include <cstdio>
  #include <cstring>
  #include <memory>
  #include <random>
  
  #include <algorithm>
  
  #include <faiss/FaissHook.h>
  #include <faiss/impl/FaissAssert.h>
  #include <faiss/utils/Heap.h>
  #include <faiss/utils/distances.h>
  #include <faiss/utils/hamming.h>
  #include <faiss/utils/utils.h>
  
  extern "C" {
  
  // general matrix multiplication
  int sgemm_(
          const char* transa,
          const char* transb,
          FINTEGER* m,
          FINTEGER* n,
          FINTEGER* k,
          const float* alpha,
          const float* a,
          FINTEGER* lda,
          const float* b,
          FINTEGER* ldb,
          float* beta,
          float* c,
          FINTEGER* ldc);
  }
  
  namespace faiss {
  
  AdditiveQuantizer::AdditiveQuantizer(
          size_t d,
          const std::vector<size_t>& nbits,
          Search_type_t search_type)
          : d(d),
            M(nbits.size()),
            nbits(nbits),
            verbose(false),
            is_trained(false),
            search_type(search_type) {
      norm_max = norm_min = NAN;
      code_size = 0;
      tot_bits = 0;
      total_codebook_size = 0;
      only_8bit = false;
      set_derived_values();
  }
  
  AdditiveQuantizer::AdditiveQuantizer()
          : AdditiveQuantizer(0, std::vector<size_t>()) {}
  
  void AdditiveQuantizer::set_derived_values() {
      tot_bits = 0;
      only_8bit = true;
      codebook_offsets.resize(M + 1, 0);
      for (int i = 0; i < M; i++) {
          int nbit = nbits[i];
          size_t k = 1 << nbit;
          codebook_offsets[i + 1] = codebook_offsets[i] + k;
          tot_bits += nbit;
          if (nbit != 0) {
              only_8bit = false;
          }
      }
      total_codebook_size = codebook_offsets[M];
      switch (search_type) {
          case ST_decompress:
          case ST_LUT_nonorm:
          case ST_norm_from_LUT:
              break; // nothing to add
          case ST_norm_float:
              tot_bits += 32;
              break;
          case ST_norm_qint8:
          case ST_norm_cqint8:
              tot_bits += 8;
              break;
          case ST_norm_qint4:
          case ST_norm_cqint4:
              tot_bits += 4;
              break;
      }
  
      // convert bits to bytes
      code_size = (tot_bits + 7) / 8;
  }
  
  namespace {
  
  // TODO
  // https://stackoverflow.com/questions/31631224/hacks-for-clamping-integer-to-0-255-and-doubles-to-0-0-1-0
  
  uint8_t encode_qint8(float x, float amin, float amax) {
      float x1 = (x - amin) / (amax - amin) * 256;
      int32_t xi = int32_t(floor(x1));
  
      return xi < 0 ? 0 : xi > 255 ? 255 : xi;
  }
  
  uint8_t encode_qint4(float x, float amin, float amax) {
      float x1 = (x - amin) / (amax - amin) * 16;
      int32_t xi = int32_t(floor(x1));
  
      return xi < 0 ? 0 : xi > 15 ? 15 : xi;
  }
  
  float decode_qint8(uint8_t i, float amin, float amax) {
      return (i + 0.5) / 256 * (amax - amin) + amin;
  }
  
  float decode_qint4(uint8_t i, float amin, float amax) {
      return (i + 0.5) / 16 * (amax - amin) + amin;
  }
  
  } // anonymous namespace
  
  uint32_t AdditiveQuantizer::encode_qcint(float x) const {
      idx_t id;
      // qnorm.assign(idx_t(1), &x, &id, idx_t(1));
      qnorm.assign(idx_t(1), &x, &id, nullptr);
      return uint32_t(id);
  }
  
  float AdditiveQuantizer::decode_qcint(uint32_t c) const {
      return qnorm.get_xb()[c];
  }
  
  void AdditiveQuantizer::pack_codes(
          size_t n,
          const int32_t* codes,
          uint8_t* packed_codes,
          int64_t ld_codes,
          const float* norms) const {
      if (ld_codes == -1) {
          ld_codes = M;
      }
      std::vector<float> norm_buf;
      if (search_type == ST_norm_float || search_type == ST_norm_qint4 ||
          search_type == ST_norm_qint8 || search_type == ST_norm_cqint8 ||
          search_type == ST_norm_cqint4) {
          if (!norms) {
              norm_buf.resize(n);
              std::vector<float> x_recons(n * d);
              decode_unpacked(codes, x_recons.data(), n, ld_codes);
              fvec_norms_L2sqr(norm_buf.data(), x_recons.data(), d, n);
              norms = norm_buf.data();
          }
      }
  #pragma omp parallel for if (n > 1000)
      for (int64_t i = 0; i < n; i++) {
          const int32_t* codes1 = codes + i * ld_codes;
          BitstringWriter bsw(packed_codes + i * code_size, code_size);
          for (int m = 0; m < M; m++) {
              bsw.write(codes1[m], nbits[m]);
          }
          switch (search_type) {
              case ST_decompress:
              case ST_LUT_nonorm:
              case ST_norm_from_LUT:
                  break;
              case ST_norm_float:
                  bsw.write(*(uint32_t*)&norms[i], 32);
                  break;
              case ST_norm_qint8: {
                  uint8_t b = encode_qint8(norms[i], norm_min, norm_max);
                  bsw.write(b, 8);
                  break;
              }
              case ST_norm_qint4: {
                  uint8_t b = encode_qint4(norms[i], norm_min, norm_max);
                  bsw.write(b, 4);
                  break;
              }
              case ST_norm_cqint8: {
                  uint32_t b = encode_qcint(norms[i]);
                  bsw.write(b, 8);
                  break;
              }
              case ST_norm_cqint4: {
                  uint32_t b = encode_qcint(norms[i]);
                  bsw.write(b, 4);
                  break;
              }
          }
      }
  }
  
  void AdditiveQuantizer::decode(const uint8_t* code, float* x, size_t n) const {
      FAISS_THROW_IF_NOT_MSG(
              is_trained, "The additive quantizer is not trained yet.");
  
      // standard additive quantizer decoding
  #pragma omp parallel for if (n > 1000)
      for (int64_t i = 0; i < n; i++) {
          BitstringReader bsr(code + i * code_size, code_size);
          float* xi = x + i * d;
          for (int m = 0; m < M; m++) {
              int idx = bsr.read(nbits[m]);
              const float* c = codebooks.data() + d * (codebook_offsets[m] + idx);
              if (m == 0) {
                  memcpy(xi, c, sizeof(*x) * d);
              } else {
                  fvec_add(d, xi, c, xi);
              }
          }
      }
  }
  
  void AdditiveQuantizer::decode_unpacked(
          const int32_t* code,
          float* x,
          size_t n,
          int64_t ld_codes) const {
      FAISS_THROW_IF_NOT_MSG(
              is_trained, "The additive quantizer is not trained yet.");
  
      if (ld_codes == -1) {
          ld_codes = M;
      }
  
      // standard additive quantizer decoding
  #pragma omp parallel for if (n > 1000)
      for (int64_t i = 0; i < n; i++) {
          const int32_t* codesi = code + i * ld_codes;
          float* xi = x + i * d;
          for (int m = 0; m < M; m++) {
              int idx = codesi[m];
              const float* c = codebooks.data() + d * (codebook_offsets[m] + idx);
              if (m == 0) {
                  memcpy(xi, c, sizeof(*x) * d);
              } else {
                  fvec_add(d, xi, c, xi);
              }
          }
      }
  }
  
  AdditiveQuantizer::~AdditiveQuantizer() {}
  
  /****************************************************************************
   * Support for fast distance computations in centroids
   ****************************************************************************/
  
  void AdditiveQuantizer::compute_centroid_norms(float* norms) const {
      size_t ntotal = (size_t)1 << tot_bits;
      // TODO: make tree of partial sums
  #pragma omp parallel
      {
          std::vector<float> tmp(d);
  #pragma omp for
          for (int64_t i = 0; i < ntotal; i++) {
              decode_64bit(i, tmp.data());
              norms[i] = fvec_norm_L2sqr(tmp.data(), d);
          }
      }
  }
  
  void AdditiveQuantizer::decode_64bit(idx_t bits, float* xi) const {
      for (int m = 0; m < M; m++) {
          idx_t idx = bits & (((size_t)1 << nbits[m]) - 1);
          bits >>= nbits[m];
          const float* c = codebooks.data() + d * (codebook_offsets[m] + idx);
          if (m == 0) {
              memcpy(xi, c, sizeof(*xi) * d);
          } else {
              fvec_add(d, xi, c, xi);
          }
      }
  }
  
  void AdditiveQuantizer::compute_LUT(size_t n, const float* xq, float* LUT)
          const {
      // in all cases, it is large matrix multiplication
  
      FINTEGER ncenti = total_codebook_size;
      FINTEGER di = d;
      FINTEGER nqi = n;
      float one = 1, zero = 0;
  
      sgemm_("Transposed",
             "Not transposed",
             &ncenti,
             &nqi,
             &di,
             &one,
             codebooks.data(),
             &di,
             xq,
             &di,
             &zero,
             LUT,
             &ncenti);
  }
  
  namespace {
  
  void compute_inner_prod_with_LUT(
          const AdditiveQuantizer& aq,
          const float* LUT,
          float* ips) {
      size_t prev_size = 1;
      for (int m = 0; m < aq.M; m++) {
          const float* LUTm = LUT + aq.codebook_offsets[m];
          int nb = aq.nbits[m];
          size_t nc = (size_t)1 << nb;
  
          if (m == 0) {
              memcpy(ips, LUT, sizeof(*ips) * nc);
          } else {
              for (int64_t i = nc - 1; i >= 0; i--) {
                  float v = LUTm[i];
                  fvec_add(prev_size, ips, v, ips + i * prev_size);
              }
          }
          prev_size *= nc;
      }
  }
  
  } // anonymous namespace
  
  void AdditiveQuantizer::knn_centroids_inner_product(
          idx_t n,
          const float* xq,
          idx_t k,
          float* distances,
          idx_t* labels) const {
      std::unique_ptr<float[]> LUT(new float[n * total_codebook_size]);
      compute_LUT(n, xq, LUT.get());
      size_t ntotal = (size_t)1 << tot_bits;
  
  #pragma omp parallel if (n > 100)
      {
          std::vector<float> dis(ntotal);
  #pragma omp for
          for (idx_t i = 0; i < n; i++) {
              const float* LUTi = LUT.get() + i * total_codebook_size;
              compute_inner_prod_with_LUT(*this, LUTi, dis.data());
              float* distances_i = distances + i * k;
              idx_t* labels_i = labels + i * k;
              minheap_heapify(k, distances_i, labels_i);
              minheap_addn(k, distances_i, labels_i, dis.data(), nullptr, ntotal);
              minheap_reorder(k, distances_i, labels_i);
          }
      }
  }
  
  void AdditiveQuantizer::knn_centroids_L2(
          idx_t n,
          const float* xq,
          idx_t k,
          float* distances,
          idx_t* labels,
          const float* norms) const {
      std::unique_ptr<float[]> LUT(new float[n * total_codebook_size]);
      compute_LUT(n, xq, LUT.get());
      std::unique_ptr<float[]> q_norms(new float[n]);
      fvec_norms_L2sqr(q_norms.get(), xq, d, n);
      size_t ntotal = (size_t)1 << tot_bits;
  
  #pragma omp parallel if (n > 100)
      {
          std::vector<float> dis(ntotal);
  #pragma omp for
          for (idx_t i = 0; i < n; i++) {
              const float* LUTi = LUT.get() + i * total_codebook_size;
              float* distances_i = distances + i * k;
              idx_t* labels_i = labels + i * k;
  
              compute_inner_prod_with_LUT(*this, LUTi, dis.data());
  
              // update distances using
              // ||x - y||^2 = ||x||^2 + ||y||^2 - 2 * <x,y>
  
              maxheap_heapify(k, distances_i, labels_i);
              for (idx_t j = 0; j < ntotal; j++) {
                  float disj = q_norms[i] + norms[j] - 2 * dis[j];
                  if (disj < distances_i[0]) {
                      heap_replace_top<CMax<float, int64_t>>(
                              k, distances_i, labels_i, disj, j);
                  }
              }
              maxheap_reorder(k, distances_i, labels_i);
          }
      }
  }
  
  /****************************************************************************
   * Support for fast distance computations in codes
   ****************************************************************************/
  
  namespace {
  
  float accumulate_IPs(
          const AdditiveQuantizer& aq,
          BitstringReader& bs,
          const uint8_t* codes,
          const float* LUT) {
      float accu = 0;
      for (int m = 0; m < aq.M; m++) {
          size_t nbit = aq.nbits[m];
          int idx = bs.read(nbit);
          accu += LUT[idx];
          LUT += (uint64_t)1 << nbit;
      }
      return accu;
  }
  
  } // anonymous namespace
  
  template <>
  float AdditiveQuantizer::
          compute_1_distance_LUT<true, AdditiveQuantizer::ST_LUT_nonorm>(
                  const uint8_t* codes,
                  const float* LUT) const {
      BitstringReader bs(codes, code_size);
      return accumulate_IPs(*this, bs, codes, LUT);
  }
  
  template <>
  float AdditiveQuantizer::
          compute_1_distance_LUT<false, AdditiveQuantizer::ST_LUT_nonorm>(
                  const uint8_t* codes,
                  const float* LUT) const {
      BitstringReader bs(codes, code_size);
      return -accumulate_IPs(*this, bs, codes, LUT);
  }
  
  template <>
  float AdditiveQuantizer::
          compute_1_distance_LUT<false, AdditiveQuantizer::ST_norm_float>(
                  const uint8_t* codes,
                  const float* LUT) const {
      BitstringReader bs(codes, code_size);
      float accu = accumulate_IPs(*this, bs, codes, LUT);
      uint32_t norm_i = bs.read(32);
      float norm2 = *(float*)&norm_i;
      return norm2 - 2 * accu;
  }
  
  template <>
  float AdditiveQuantizer::
          compute_1_distance_LUT<false, AdditiveQuantizer::ST_norm_cqint8>(
                  const uint8_t* codes,
                  const float* LUT) const {
      BitstringReader bs(codes, code_size);
      float accu = accumulate_IPs(*this, bs, codes, LUT);
      uint32_t norm_i = bs.read(8);
      float norm2 = decode_qcint(norm_i);
      return norm2 - 2 * accu;
  }
  
  template <>
  float AdditiveQuantizer::
          compute_1_distance_LUT<false, AdditiveQuantizer::ST_norm_cqint4>(
                  const uint8_t* codes,
                  const float* LUT) const {
      BitstringReader bs(codes, code_size);
      float accu = accumulate_IPs(*this, bs, codes, LUT);
      uint32_t norm_i = bs.read(4);
      float norm2 = decode_qcint(norm_i);
      return norm2 - 2 * accu;
  }
  
  template <>
  float AdditiveQuantizer::
          compute_1_distance_LUT<false, AdditiveQuantizer::ST_norm_qint8>(
                  const uint8_t* codes,
                  const float* LUT) const {
      BitstringReader bs(codes, code_size);
      float accu = accumulate_IPs(*this, bs, codes, LUT);
      uint32_t norm_i = bs.read(8);
      float norm2 = decode_qint8(norm_i, norm_min, norm_max);
      return norm2 - 2 * accu;
  }
  
  template <>
  float AdditiveQuantizer::
          compute_1_distance_LUT<false, AdditiveQuantizer::ST_norm_qint4>(
                  const uint8_t* codes,
                  const float* LUT) const {
      BitstringReader bs(codes, code_size);
      float accu = accumulate_IPs(*this, bs, codes, LUT);
      uint32_t norm_i = bs.read(4);
      float norm2 = decode_qint4(norm_i, norm_min, norm_max);
      return norm2 - 2 * accu;
  }
  
  } // namespace faiss
  