gvarHash.hpp

Go to the documentation of this file.

/*
 * Copyright (c) 2021 Anthony J. Greenberg
 *
 * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
 * IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
 * THE POSSIBILITY OF SUCH DAMAGE.
 */
 
 
#pragma once
 
#include <cstddef>
#include <fstream>
#include <vector>
#include <utility>        // for std::pair
#include <string>
#include <thread>
 
#include "similarityMatrix.hpp"
 
namespace BayesicSpace {
    struct LocationWithLength;
    struct CountAndSize;
    struct IndividualAndSketchCounts;
    struct BedDataStats;
    struct InOutFileNames;
    struct SparsityParameters;
    struct HashGroup;
    struct HashGroupItPairCount;
    class GenoTableBin;
    class GenoTableHash;
 
    struct LocationWithLength {
        size_t start;
        size_t length;
    };
 
    struct CountAndSize {
        size_t count;
        size_t size;
    };
 
    struct IndividualAndSketchCounts {
        uint32_t nIndividuals;
        uint16_t kSketches;
    };
 
    struct BedDataStats {
        size_t firstLocusIdx;
        size_t nLociPerThread;
        size_t nBytesPerLocus;
        size_t nBytesToRead;
        size_t nLociToRead;
        size_t nMemChunks;
    };
    struct InOutFileNames {
        std::string inputFileName;
        std::string outputFileName;
    };
    struct SparsityParameters{
        size_t nRowsPerBand;
        float similarityCutOff;
    };
    struct HashGroup {
        uint64_t cumulativeNpairs;
        std::vector<uint32_t> locusIndexes;
    };
 
    struct HashGroupItPairCount {
        size_t pairCount{0};
        std::vector<HashGroup>::const_iterator hgIterator;
    };
 
    class GenoTableBin {
    public:
        GenoTableBin() : nIndividuals_{0}, nLoci_{0}, binLocusSize_{0}, nThreads_{1} {};
        GenoTableBin(const std::string &inputFileName, const uint32_t &nIndividuals, std::string logFileName) : 
                        GenoTableBin( inputFileName, nIndividuals, std::move(logFileName), std::thread::hardware_concurrency() ) {};
        GenoTableBin(const std::string &inputFileName, const uint32_t &nIndividuals, std::string logFileName, const size_t &nThreads);
        GenoTableBin(const std::vector<int> &maCounts, const uint32_t &nIndividuals, std::string logFileName) :
                        GenoTableBin( maCounts, nIndividuals, std::move(logFileName), std::thread::hardware_concurrency() ) {};
        GenoTableBin(const std::vector<int> &maCounts, const uint32_t &nIndividuals, std::string logFileName, const size_t &nThreads);
 
        GenoTableBin(const GenoTableBin &toCopy) = delete;
        GenoTableBin& operator=(const GenoTableBin &toCopy) = delete;
        GenoTableBin(GenoTableBin &&toMove) noexcept = default;
        GenoTableBin& operator=(GenoTableBin &&toMove) noexcept = default;
        ~GenoTableBin() = default;
 
        void saveGenoBinary(const std::string &outFileName) const;
        void allJaccardLD( const InOutFileNames &bimAndLDnames, const size_t &suggestNchunks = static_cast<size_t>(1) ) const;
        void saveLogFile() const;
    private:
        std::vector<uint8_t> binGenotypes_;
        uint32_t nIndividuals_;
        uint32_t nLoci_;
        size_t binLocusSize_;
        size_t nThreads_;
        static const size_t nMagicBytes_;
        static const uint8_t oneBit_;
        static const uint8_t byteSize_;
        static const uint8_t bedGenoPerByte_;
        static const uint8_t llWordSize_;
        mutable std::string logMessages_;
        std::string logFileName_;
        void bed2binBlk_(const std::vector<char> &bedData, const std::pair<size_t, size_t> &bedLocusIndRange, const LocationWithLength &locusSpan);
        size_t bed2binThreaded_(const std::vector<char> &bedData, const std::vector< std::pair<size_t, size_t> > &threadRanges, const LocationWithLength &locusSpan);
        size_t bed2bin_(const BedDataStats &locusGroupStats, std::fstream &bedStream);
        void mac2binBlk_(const std::vector<int> &macData, const std::pair<size_t, size_t> &locusIndRange);
        [[gnu::warn_unused_result]] SimilarityMatrix jaccardBlock_(const std::pair<RowColIdx, RowColIdx> &blockRange) const;
        [[gnu::warn_unused_result]] SimilarityMatrix jaccardThreaded_(const std::vector< std::pair<RowColIdx, RowColIdx> > &indexPairs) const;
        [[gnu::warn_unused_result]] JaccardPair makeJaccardPair_(const RowColIdx &rowColumn) const;
    };
    class GenoTableHash {
    public:
        GenoTableHash() : nIndividuals_{0}, kSketches_{0}, sketchSize_{0}, nLoci_{0}, locusSize_{0}, nFullWordBytes_{0}, nThreads_{1}, emptyBinIdxSeed_{0} {};
        GenoTableHash(const std::string &inputFileName, const IndividualAndSketchCounts &indivSketchCounts, const size_t &nThreads, std::string logFileName);
        GenoTableHash(const std::string &inputFileName, const IndividualAndSketchCounts &indivSketchCounts, std::string logFileName) :
                            GenoTableHash( inputFileName, indivSketchCounts, std::thread::hardware_concurrency(), std::move(logFileName) ) {};
        GenoTableHash(const std::vector<int> &maCounts, const IndividualAndSketchCounts &indivSketchCounts, const size_t &nThreads, std::string logFileName);
        GenoTableHash(const std::vector<int> &maCounts, const IndividualAndSketchCounts &indivSketchCounts, std::string logFileName) :
                GenoTableHash( maCounts, indivSketchCounts, std::thread::hardware_concurrency(), std::move(logFileName) ) {};
 
        GenoTableHash(const GenoTableHash &toCopy) = delete;
        GenoTableHash& operator=(const GenoTableHash &toCopy) = delete;
        GenoTableHash(GenoTableHash &&toMove) noexcept = default;
        GenoTableHash& operator=(GenoTableHash &&toMove) noexcept = default;
        ~GenoTableHash() = default;
 
        void allHashLD( const float &similarityCutOff, const InOutFileNames &bimAndLDnames, const size_t &suggestNchunks = static_cast<size_t>(1) ) const;
        [[gnu::warn_unused_result]] std::vector<HashGroup> makeLDgroups(const size_t &nRowsPerBand) const;
        void makeLDgroups(const size_t &nRowsPerBand, const InOutFileNames &bimAndGroupNames) const;
        void ldInGroups(const SparsityParameters &sparsityValues, const InOutFileNames &bimAndLDnames, const size_t &suggestNchunks = static_cast<size_t>(1) ) const;
        void saveLogFile() const;
    private:
        std::vector<uint16_t> sketches_;
        uint32_t nIndividuals_;
        uint16_t kSketches_;
        uint32_t sketchSize_;
        uint32_t nLoci_;
        uint32_t locusSize_;
        size_t nFullWordBytes_;
        size_t nThreads_;
        uint64_t emptyBinIdxSeed_;
        mutable std::string logMessages_;
        std::string logFileName_;
        static const size_t nMagicBytes_;
        static const uint8_t oneBit_;
        static const uint8_t byteSize_;
        static const uint8_t bedGenoPerByte_;
        static const uint8_t llWordSize_;
        static const uint32_t roundMask_;
        static const uint64_t allBitsSet_;
        static const size_t wordSizeInBits_;
        static const uint16_t emptyBinToken_;
        void permuteBits_(const std::vector<size_t> &permutationIdx, std::vector<uint8_t> &binLocus) const;
        void locusOPH_(const size_t &locusInd, const std::vector<size_t> &permutation, std::vector<uint8_t> &binLocus);
        void bed2ophBlk_(const std::vector<char> &bedData, const std::pair<size_t, size_t> &bedLocusIndRange, const LocationWithLength &bedLocusSpan,
                                    const std::vector<size_t> &permutation, const std::vector< std::pair<size_t, size_t> > &padIndiv);
        size_t bed2ophThreaded_(const std::vector<char> &bedData, const std::vector< std::pair<size_t, size_t> > &threadRanges, const LocationWithLength &bedLocusSpan,
                            const std::vector<size_t> &permutation, const std::vector< std::pair<size_t, size_t> > &padIndiv);
        size_t bed2oph_(const BedDataStats &locusGroupStats, std::fstream &bedStream, const std::vector<size_t> &permutation,
                            const std::vector< std::pair<size_t, size_t> > &padIndiv);
        void mac2ophBlk_(const std::vector<int> &macData, const std::pair<size_t, size_t> &blockRange,
                const std::vector<size_t> &permutation, const std::vector< std::pair<size_t, size_t> > &padIndiv);
        size_t mac2ophThreaded_(const std::vector<int> &macData, const std::vector< std::pair<size_t, size_t> > &threadRanges, const LocationWithLength &locusBlock,
                            const std::vector<size_t> &permutation, const std::vector< std::pair<size_t, size_t> > &padIndiv);
        [[gnu::warn_unused_result]] SimilarityMatrix hashJacBlock_(const std::pair<RowColIdx, RowColIdx> &blockRange, const std::vector<uint32_t> &locusIndexes, const float &similarityCutOff) const;
        [[gnu::warn_unused_result]] SimilarityMatrix hashJacBlock_(const std::pair<HashGroupItPairCount, HashGroupItPairCount> &blockRange, const float &similarityCutOff) const;
        [[gnu::warn_unused_result]] SimilarityMatrix hashJacThreaded_(const std::vector< std::pair<RowColIdx, RowColIdx> > &indexPairs,
                const std::vector<uint32_t> &locusIndexes, const float &similarityCutOff) const;
        [[gnu::warn_unused_result]] SimilarityMatrix hashJacThreaded_(const std::vector< std::pair<HashGroupItPairCount, HashGroupItPairCount> > &blockRanges, const float &similarityCutOff) const;
        [[gnu::warn_unused_result]] JaccardPair makeJaccardPair_(const RowColIdx &rowColumn) const noexcept;
    };
}