pfor2008.h

/**
 * This code is released under the
 * Apache License Version 2.0 http://www.apache.org/licenses/.
 *
 * (c) Daniel Lemire, http://lemire.me/en/
 */

#ifndef PFOR2008_H_
#define PFOR2008_H_

#include "common.h"
#include "codecs.h"
#include "bitpacking.h"
#include "util.h"
#include <iostream>

namespace FastPForLib {

/**
 * Follows:
 *
 *  Zhang J, Long X, Suel T. Performance of compressed inverted list caching in search engines.
 *  Proceeedings of 17th Conference on the World Wide Web, Beijing, China, Huai J, Chen R,
 *  Hon H-W, Liu Y, Ma W-Y, Tomkins A, Zhang X (eds.). ACM Press: New York, April 2008; 387�396.
 *
 * It is very similar to PFor except that exceptions can be stored using 8, 16 or 32 bits.
 * So I just copied PFor and made the required changes.
 *
 * In a multithreaded context, you may need one PFor2008 per thread.
 *
 */
class PFor2008: public IntegerCODEC {
public:
    enum {
        BlockSizeInUnitsOfPackSize = 4,
        PACKSIZE = 32,
        BlockSize = BlockSizeInUnitsOfPackSize * PACKSIZE,
        blocksizeinbits = 7//constexprbits(BlockSize)
    };
    // these are reusable buffers
    std::vector<uint32_t> codedcopy;
    std::vector<uint32_t> miss;
    typedef uint32_t DATATYPE;// this is so that our code looks more like the original paper

    PFor2008() :
        codedcopy(BlockSize), miss(BlockSize) {
    }
    // for delta coding, we don't use a base.
    static uint32_t determineBestBase(const DATATYPE * in, size_t size,
            const uint32_t maxb) {
        if (size == 0)
            return 0;
        const size_t defaultsamplesize = 64 * 1024;
        // the original paper describes sorting
        // a sample, but this only makes sense if you
        // are coding a frame of reference.
        size_t samplesize = size > defaultsamplesize ? defaultsamplesize : size;
        std::vector<uint32_t> freqs(33);
        // we choose the sample to be consecutive
        uint32_t rstart = size > samplesize ? (rand() % (static_cast<uint32_t>(size - samplesize)))
                : 0U;
        for (uint32_t k = rstart; k < rstart + samplesize; ++k) {
            freqs[asmbits(in[k])]++;
        }
        uint32_t bestb = maxb;
        uint32_t numberofexceptions = 0;
        double Erate = 0;
        double bestcost = maxb;
        for (uint32_t b = bestb - 1; b < maxb; --b) {
            numberofexceptions += freqs[b + 1];
            Erate = static_cast<double>(numberofexceptions) / static_cast<double>(samplesize);
            /**
             * though this is not explicit in the original paper, you
             * need to somehow compensate for compulsory exceptions
             * when the chosen number of bits is small.
             *
             * We use their formula (3.1.5) to estimate actual number
             * of total exceptions, including compulsory exceptions.
             */
            if (numberofexceptions > 0) {
                double altErate = (Erate * 128 - 1) / (Erate * (1U << b));
                if (altErate > Erate)
                    Erate = altErate;
            }
            const double thiscost = b + Erate * maxb;
            if (thiscost <= bestcost) {
                bestcost = thiscost;
                bestb = b;
            }
        }
        return bestb;
    }

    // returns location of first exception or BlockSize if there is none
    uint32_t compressblockPFOR(const DATATYPE * __restrict__ in,
            uint32_t * __restrict__ outputbegin, const uint32_t b,
            DATATYPE * __restrict__ & exceptions) {
        if (b == 32) {
            for (size_t k = 0; k < BlockSize; ++k)
                *(outputbegin++) = *(in++);
            return BlockSize;
        }
        size_t exceptcounter = 0;
        const uint32_t maxgap = 1U << b;
        {
            std::vector<uint32_t>::iterator cci = codedcopy.begin();
            for (uint32_t k = 0; k < BlockSize; ++k, ++cci) {
                miss[exceptcounter] = k;
                exceptcounter += (in[k] >= maxgap);
            }

        }
        if (exceptcounter == 0) {
            packblock(in, outputbegin, b);
            return BlockSize;
        }
        codedcopy.assign(in, in + BlockSize);
        uint32_t firstexcept = miss[0];
        uint32_t prev = 0;
        *(exceptions++) = codedcopy[firstexcept];
        prev = firstexcept;
        if (maxgap < BlockSize) {
            for (size_t i = 1; i < exceptcounter; ++i) {
            	uint32_t cur = miss[i];
                // they don't include this part, but it is required:
                while (cur > maxgap + prev) {
                    // compulsory exception
                    uint32_t compulcur = prev + maxgap;
                    *(exceptions++) = codedcopy[compulcur];
                    codedcopy[prev] = maxgap - 1;
                    prev = compulcur;
                }
                *(exceptions++) = codedcopy[cur];
                codedcopy[prev] = cur - prev - 1;
                prev = cur;
            }
        } else {
            for (uint32_t i = 1; i < exceptcounter; ++i) {
                uint32_t cur = miss[i];
                *(exceptions++) = codedcopy[cur];
                codedcopy[prev] = cur - prev - 1;
                prev = cur;
            }
        }
        packblock(&codedcopy[0], outputbegin, b);
        return firstexcept;
    }

    void packblock(const uint32_t * source, uint32_t * out, const uint32_t bit) {
        for (uint32_t j = 0; j != BlockSize; j += PACKSIZE) {
            fastpack(source + j, out, bit);
            out += bit;
        }
    }

    void unpackblock(const uint32_t * source, uint32_t * out,
            const uint32_t bit) {
        for (uint32_t j = 0; j != BlockSize; j += PACKSIZE) {
            fastunpack(source, out + j, bit);
            source += bit;
        }
    }

    void encodeArray(const uint32_t *in, const size_t len, uint32_t *out,
            size_t &nvalue) {
        *out++ = static_cast<uint32_t>(len);
#ifndef NDEBUG
        const uint32_t * const finalin(in + len);
#endif
        const uint32_t maxsize = (1U << (32 - blocksizeinbits - 1));
        size_t totalnvalue(1);
        //for (size_t i = 0; i < len; i += maxsize)
        for (size_t j = 0; j < (len+maxsize-1U)/maxsize; ++j) {
            size_t i = j << (32 - blocksizeinbits - 1);
            size_t l = maxsize;
            if (i + maxsize > len) {
                l = len - i;
                assert(l <= maxsize);
            }
            size_t thisnvalue = nvalue - totalnvalue;
            assert(in + i + l <= finalin);
            __encodeArray(&in[i], l, out, thisnvalue);
            totalnvalue += thisnvalue;
            assert(totalnvalue <= nvalue);
            out += thisnvalue;
        }
        nvalue = totalnvalue;
    }
    const uint32_t * decodeArray(const uint32_t *in, const size_t len,
            uint32_t *out, size_t &nvalue) {
        nvalue = *in++;
        if (nvalue == 0)
            return in;
#ifndef NDEBUG
        const uint32_t * const initin = in;
#endif
        const uint32_t * const finalin = in + len;
        size_t totalnvalue(0);
        while (totalnvalue < nvalue) {
            size_t thisnvalue = nvalue - totalnvalue;
#ifndef NDEBUG
            const uint32_t * const befin(in);
#endif
            assert(finalin <= len + in);
            const uint32_t maxb = *in++;
            if (maxb == 32) {
                in
                        = __decodeArray<uint32_t> (in, finalin - in, out,
                                thisnvalue);
            } else if (maxb == 16) {
                in
                        = __decodeArray<uint16_t> (in, finalin - in, out,
                                thisnvalue);
            } else if (maxb == 8) {
                in = __decodeArray<uint8_t> (in, finalin - in, out, thisnvalue);
            } else {
                throw std::logic_error("corrupted?");
            }
            assert(in > befin);
            assert(in <= finalin);
            out += thisnvalue;
            totalnvalue += thisnvalue;
            assert(totalnvalue <= nvalue);
        }
        assert(in <= len + initin);
        assert(in <= finalin);
        nvalue = totalnvalue;
        return in;

    }

    uint32_t howmanybits(const uint32_t *in, const size_t len) {
        uint32_t accumulator = 0;
        for (uint32_t k = 0; k < len; ++k) {
            accumulator |= in[k];
        }
        if (accumulator >= (1U << 16))
            return 32;
        if (accumulator >= (1U << 8))
            return 16;
        return 8;
    }

    void __encodeArray(const uint32_t *in, const size_t len, uint32_t *out,
            size_t &nvalue) {
        const uint32_t maxb = howmanybits(in, len);
        checkifdivisibleby(len, BlockSize);
        const uint32_t * const initout(out);
        std::vector<DATATYPE> exceptions(len);
        DATATYPE * __restrict__ i = &exceptions[0];
        const uint32_t b = determineBestBase(in, len, maxb);
        *out++ = maxb;
        *out++ = static_cast<uint32_t>(len);
        *out++ = b;

        for (size_t k = 0; k < len / BlockSize; ++k) {
            uint32_t * const headerout(out);
            ++out;
            uint32_t firstexcept = compressblockPFOR(in, out, b, i);
            out += (BlockSize * b) / 32;
            in += BlockSize;
            const uint32_t bitsforfirstexcept = blocksizeinbits;
            const uint32_t firstexceptmask = (1U << blocksizeinbits) - 1;
            const uint32_t exceptindex = static_cast<uint32_t>(i - &exceptions[0]);
            *headerout = (firstexcept & firstexceptmask) | (exceptindex
                    << bitsforfirstexcept);
        }
        const size_t howmanyexcept = i - &exceptions[0];
        if (maxb == 32) {
            for (uint32_t t = 0; t < howmanyexcept; ++t)
                *out++ = exceptions[t];
        } else if (maxb == 16) {
            uint16_t * out16 = reinterpret_cast<uint16_t *> (out);
            for (uint32_t t = 0; t < howmanyexcept; ++t)
                *out16++ = static_cast<uint16_t>(exceptions[t]);
            out
                    = reinterpret_cast<uint32_t *> ((reinterpret_cast<uintptr_t> (out16)
                            + 3) & ~3);
        } else if (maxb == 8) {
            uint8_t * out8 = reinterpret_cast<uint8_t *> (out);
            for (uint32_t t = 0; t < howmanyexcept; ++t)
                *out8++ = static_cast<uint8_t>(exceptions[t]);
            out
                    = reinterpret_cast<uint32_t *> ((reinterpret_cast<uintptr_t> (out8)
                            + 3) & ~3);
        } else {
            throw std::logic_error("should not happen");
        }

        nvalue = out - initout;
    }
    template<class EXCEPTTYPE>
#ifndef NDEBUG
    const uint32_t * __decodeArray(const uint32_t *in, const size_t len,
#else
    const uint32_t * __decodeArray(const uint32_t *in, const size_t    ,
#endif
            uint32_t *out, size_t &nvalue) {
#ifndef NDEBUG
        const uint32_t * const initin(in);
#endif
        nvalue = *in++;
        checkifdivisibleby(nvalue, BlockSize);
        const uint32_t b = *in++;
        const EXCEPTTYPE * __restrict__ except =
                reinterpret_cast<const EXCEPTTYPE *> (in + nvalue * b / 32
                        + nvalue / BlockSize);
        const uint32_t bitsforfirstexcept = blocksizeinbits;
        const uint32_t firstexceptmask = (1U << blocksizeinbits) - 1;
        const EXCEPTTYPE * endexceptpointer = except;
        const EXCEPTTYPE * const initexcept(except);
        for (size_t k = 0; k < nvalue / BlockSize; ++k) {
            const uint32_t * const headerin(in);
            ++in;
            const uint32_t firstexcept = *headerin & firstexceptmask;
            const uint32_t exceptindex = *headerin >> bitsforfirstexcept;
            endexceptpointer = initexcept + exceptindex;
            uncompressblockPFOR(in, out, b, except, endexceptpointer,
                    firstexcept);
            in += (BlockSize * b) / 32;
            out += BlockSize;
        }
        assert(initin + len >= in);
        assert(
                initin + len
                        >= reinterpret_cast<const uint32_t *> (endexceptpointer));
        return reinterpret_cast<const uint32_t *> ((reinterpret_cast<uintptr_t> (endexceptpointer)
                + 3) & ~3);
    }

    template<class EXCEPTTYPE>
    void uncompressblockPFOR(
            const uint32_t * __restrict__ inputbegin,// points to the first packed word
            DATATYPE * __restrict__ outputbegin,
            const uint32_t b,
            const EXCEPTTYPE *__restrict__ & i, // i points to value of the first exception
            const EXCEPTTYPE * __restrict__ end_exception,
            size_t next_exception // points to the position of the first exception
    ) {
        unpackblock(inputbegin, reinterpret_cast<uint32_t *> (outputbegin), b); /* bit-unpack the values */
        for (size_t cur = next_exception; i != end_exception; cur
                = next_exception) {
            next_exception = cur + static_cast<size_t> (outputbegin[cur]) + 1;
            outputbegin[cur] = *(i++);
        }
    }

    virtual std::string name() const {
        std::ostringstream convert;
        convert << "PFor2008";
        return convert.str();
    }

};

} // namespace FastPFor

#endif /* PFOR2008_H_ */