strat.h

Go to the documentation of this file.

/******************************************************************************
 *
 * Project: sgs
 * Purpose: C++ implementation of stratified sampling
 * Author: Joseph Meyer
 * Date: September, 2025
 *
 ******************************************************************************/

 
#include <iostream>
#include <random>
 
#include "utils/access.h"
#include "utils/existing.h"
#include "utils/helper.h"
#include "utils/raster.h"
#include "utils/vector.h"
 
#include <xoshiro.h>
 
namespace sgs {
namespace strat {

std::vector<int64_t>
calculateAllocation(
    int64_t numSamples,
    std::string allocation, 
    std::vector<int64_t> strataCounts, 
    std::vector<double> weights,
    int64_t numPixels)
{
    std::vector<int64_t> retval;
    int64_t remainder = numSamples;
    int64_t numStrata = strataCounts.size();
    if (allocation == "prop") {
        //allocate the samples per stratum according to stratum size
        int64_t pixelsPerSample = numPixels / numSamples;
 
        //add 1 if pixelsPerSample was truncated down, to avoid adding too many samples
        pixelsPerSample += static_cast<int64_t>(pixelsPerSample * numSamples < numPixels);
 
        for (int64_t i = 0; i < numStrata; i++) {
            int64_t count = strataCounts[i] / pixelsPerSample;
            retval.push_back(count);
            remainder -= count;
        }
    }
    else if (allocation == "equal") {
        //determine the count of samples per strata
        int64_t strataSampleCount = numSamples / numStrata;
        
        for (int64_t i = 0; i < numStrata; i++) {
            retval.push_back(strataSampleCount);
            remainder -= strataSampleCount;
        }
    }
    else if (allocation == "manual" || allocation == "optim") { 
        //allocate samples accordign to weights. 
        //Optim allocation calculates these weights whereas when using manual the weights are given by the user
        for (int64_t i = 0; i < numStrata; i++) {
            int64_t count = static_cast<int64_t>(static_cast<double>(numSamples) * weights[i]);
            retval.push_back(count);
            remainder -= count;
        }
    }
    else { 
        throw std::runtime_error("allocation method must be one of 'prop', 'equal', 'manual', or 'optim'.");
    }
 
    //redistribute remainder pixels among strata, and check strata sizes
    for (int64_t i = numStrata; i > 0; i--) {
        int64_t extra = remainder / i;
        retval[i - 1] += extra;
        remainder -= extra;
 
        if (retval[i - 1] > strataCounts[i - 1]) {
            std::cout << "warning: strata " << i - 1 << " does not have enough pixels for the full " << retval[i - 1] << " samples it should recieve. There will be less than " << numSamples << " final samples." << std::endl;
            retval[i - 1] = strataCounts[i - 1];
        }
    }
 
    return retval;
}

struct OptimAllocationDataManager {
    helper::RasterBandMetaData band;
    std::vector<helper::Variance> variances;
    bool used = false;

    OptimAllocationDataManager(raster::GDALRasterWrapper *p_raster, int bandNum, std::string allocation) {
        if (!p_raster || allocation != "optim") {
            return;
        }
 
        this->band.p_band = p_raster->getRasterBand(bandNum);
        this->band.type = p_raster->getRasterBandType(bandNum);
        this->band.size = p_raster->getRasterBandTypeSize(bandNum);
        this->band.p_buffer = nullptr; //allocated later
        this->band.nan = this->band.p_band->GetNoDataValue();
        this->band.p_band->GetBlockSize(&this->band.xBlockSize, &this->band.yBlockSize);
        this->used = true;
    }
    
    ~OptimAllocationDataManager() {
        if (this->used) {
            VSIFree(this->band.p_buffer);
        }
    }

    inline void
    init(int numStrata, int xBlockSize, int yBlockSize) {
        this->variances.resize(numStrata);
        this->band.p_buffer = VSIMalloc3(xBlockSize, yBlockSize, band.size);
    }

    inline void
    readNewBlock(int xBlockSize, int yBlockSize, int xBlock, int yBlock, int xValid, int yValid) {
        helper::rasterBandIO(this->band, this->band.p_buffer, xBlockSize, yBlockSize, xBlock, yBlock, xValid, yValid, true, false);
    }

    inline void
    update(int index, int strata) {
        double val = helper::getPixelValueDependingOnType<double>(this->band.type, this->band.p_buffer, index);
        variances[strata].update(val); 
    }

    inline std::vector<double>
    getAllocationPercentages(void) {
        std::vector<double> retval(variances.size());
 
        double total = 0;
        for (size_t i = 0; i < variances.size(); i++) {
            double stdev = variances[i].getStdev();
            double count = static_cast<double>(variances[i].getCount());
            double product = count == 0 ? 0 : stdev * count; //if count == 0 stdev will be nan, this stops the nan from spreading
            retval[i] = product;
            total += product;
        }
 
        for (size_t i = 0; i < variances.size(); i++) {
            retval[i] = retval[i] / total;
        }
 
        return retval;
    }
};

class IndexStorageVectors {
private:
    std::vector<int64_t> strataCounts;
 
    std::vector<int64_t> indexCountPerStrata;
    std::vector<std::vector<helper::Index>> indexesPerStrata;
    
    std::vector<int64_t> firstXIndexCountPerStrata;
    std::vector<std::vector<helper::Index>> firstXIndexesPerStrata;
 
    int64_t numStrata;
    int64_t x;
 
public:
    IndexStorageVectors(int64_t numStrata, int64_t x) {
        this->numStrata = numStrata;
        this->x = x;
 
        this->strataCounts.resize(numStrata);
        
        this->indexCountPerStrata.resize(numStrata, 0);
        this->indexesPerStrata.resize(numStrata);
        
        this->firstXIndexCountPerStrata.resize(numStrata, 0);
        this->firstXIndexesPerStrata.resize(numStrata);
        for (int64_t i = 0; i < numStrata; i++) {
            this->firstXIndexesPerStrata[i].resize(x);
        }
    }

    inline void
    updateStrataCounts(int strata) {
        this->strataCounts[strata]++;
    }

    inline void
    updateIndexesVector(int strata, helper::Index& index) {
        this->indexesPerStrata[strata].push_back(index);
        this->indexCountPerStrata[strata]++;
    }

    inline void
    updateFirstXIndexesVector(int strata, helper::Index& index) {
        int64_t i = firstXIndexCountPerStrata[strata];
        if (i < this->x) {
            firstXIndexesPerStrata[strata][i] = index;
            firstXIndexCountPerStrata[strata]++;    
        }
        else if (i == this->x) {
            std::vector<helper::Index>().swap(firstXIndexesPerStrata[strata]);
            firstXIndexCountPerStrata[strata]++;
        }
    }

    inline std::vector<std::vector<helper::Index> *>
    getStrataIndexVectors(std::vector<int64_t> existing, std::vector<int64_t> strataSampleCounts, xso::xoshiro_4x64_plus& rng) {
        std::vector<std::vector<helper::Index> *> retval(numStrata);
 
        for (int64_t i = 0; i < numStrata; i++) {
            int64_t existingSamples = existing[i];
            int64_t desiredSamples = strataSampleCounts[i];
            int64_t remainingSamples = desiredSamples - existingSamples;
 
            int64_t probIndexesCount = indexCountPerStrata[i];
            int64_t firstXIndexesCount = firstXIndexCountPerStrata[i];
 
            if (probIndexesCount >= remainingSamples || firstXIndexesCount > x) {
                auto begin = this->indexesPerStrata[i].begin();
                auto end = this->indexesPerStrata[i].end();
                std::shuffle(begin, end, rng);
                retval[i] = &this->indexesPerStrata[i];
            }
            else {
                this->firstXIndexesPerStrata[i].resize(firstXIndexesCount);
                auto begin = this->firstXIndexesPerStrata[i].begin();
                auto end = this->firstXIndexesPerStrata[i].end();
                std::shuffle(begin, end, rng);
                retval[i] = &this->firstXIndexesPerStrata[i];
            }
        }
 
        return retval;
    }

    inline std::vector<int64_t>
    getStrataCounts(void) {
        return this->strataCounts;
    }

    inline int64_t
    getNumDataPixels(void) {
        int64_t retval = 0;
 
        for (const int64_t& count : this->strataCounts) {
            retval += count;
        }
 
        return retval;
    }
};

template <typename T>
std::vector<int64_t>
processBlocksStratRandom(
    int numSamples,
    int numStrata,
    helper::RasterBandMetaData& band,
    access::Access& access,
    existing::Existing& existing,
    IndexStorageVectors& indices,
    std::vector<std::vector<OGRPoint>>& existingSamples,
    uint64_t multiplier,
    xso::xoshiro_4x64_plus& rng,
    std::string allocation,
    OptimAllocationDataManager& optim,
    std::vector<double> weights,
    int width,
    int height) 
{
    T nanInt = static_cast<T>(band.nan);
    int xBlockSize = band.xBlockSize;
        int yBlockSize = band.yBlockSize;
 
    int xBlocks = (width + xBlockSize - 1) / xBlockSize;
    int yBlocks = (height + yBlockSize - 1) / yBlockSize;
 
    band.p_buffer = VSIMalloc3(xBlockSize, yBlockSize, band.size);
    T *p_buffer = reinterpret_cast<T *>(band.p_buffer);
    int8_t *p_access = nullptr;
    if (access.used) {
        access.band.p_buffer = VSIMalloc3(xBlockSize, yBlockSize, access.band.size);
        p_access = reinterpret_cast<int8_t *>(access.band.p_buffer);
    }
 
    helper::RandValController rand(band.xBlockSize, band.yBlockSize, multiplier, &rng);
 
    if (optim.used) {
        optim.init(numStrata, xBlockSize, yBlockSize);
    }
 
    for (int yBlock = 0; yBlock < yBlocks; yBlock++) {
        for (int xBlock = 0; xBlock < xBlocks; xBlock++) {
            int xValid, yValid;
    
            //read block
            band.p_band->GetActualBlockSize(xBlock, yBlock, &xValid, &yValid);
            helper::rasterBandIO(band, band.p_buffer, xBlockSize, yBlockSize, xBlock, yBlock, xValid, yValid, true, false);
 
            //read access block
            if (access.used) {
                helper::rasterBandIO(access.band, access.band.p_buffer, xBlockSize, yBlockSize, xBlock, yBlock, xValid, yValid, true, false);
            }
 
            //read mraster block for optim
            if (optim.used) {
                optim.readNewBlock(xBlockSize, yBlockSize, xBlock, yBlock, xValid, yValid);
            }
 
            //calculate rand vals
            rand.calculateRandValues();
 
            //iterate through block and update vectors
            for (int y = 0; y < yValid; y++) {
                int blockIndex = y * xBlockSize;
                for (int x = 0; x < xValid; x++) {
                    T val = p_buffer[blockIndex];
                    helper::Index index = {x + xBlock * xBlockSize, y + yBlock * yBlockSize};
 
                    bool isNan = val == nanInt;
                    bool accessible = !access.used || p_access[blockIndex] != 1;
                    bool alreadySampled = existing.used && existing.containsIndex(index.x, index.y);
 
                    //check nan
                    if (isNan) {
                        blockIndex++;
                        continue;
                    }
 
                    if (val >= numStrata) {
                        throw std::runtime_error("the num_strata indicated for the strat raster band is less than or equal to one of the value sin that band.");
                    }
 
                    if (val < 0) {
                        std::string errmsg = "a negative value of " + std::to_string(val) + " was found in the strat raster, and has not been marked as a nodata value.";
                        throw std::runtime_error(errmsg);
                    }
 
                    //update optim allocation variance calculations
                    if (optim.used) {
                        optim.update(blockIndex, val);
                    }
 
                    //udpate strata counts
                    indices.updateStrataCounts(val);
                    
                    //update existing sampled strata
                    if (alreadySampled) {
                        existingSamples[val].push_back(existing.getPoint(index.x, index.y));
                    }
 
                    //add val to stored indices
                    if (accessible && !alreadySampled) {
                        indices.updateFirstXIndexesVector(val, index);
 
                        if (rand.next()) {
                            indices.updateIndexesVector(val, index);
                        }                   
                    }
 
                    //increment block index
                    blockIndex++;
                }
            }
        }
    }
 
    VSIFree(band.p_buffer);
    if (access.used) {
        VSIFree(access.band.p_buffer);
    }   
 
    if (optim.used) {
        weights = optim.getAllocationPercentages();
    }
 
    std::vector<int64_t> strataSampleCounts = calculateAllocation(
        numSamples,
        allocation,
        indices.getStrataCounts(),
        weights,
        indices.getNumDataPixels()
    );
 
    return strataSampleCounts;
}

struct FocalWindow {
    int wrow, wcol;
    int width;
    int vpad, hpad;
    std::vector<bool> m;
    std::vector<bool> valid;

    FocalWindow(int wrow, int wcol, int width) {
        this->wrow = wrow;
        this->vpad = wrow / 2;
        this->wcol = wcol;
        this->hpad = wcol / 2;
        this->width = width;
        this->m.resize(wrow * width, false);
        this->valid.resize(wrow * width, false);
    }

    inline void
    reset(int row) {
        int start = (row % wrow) * this->width;
        int end = start + this->width;
        for (int i = start; i < end; i++) {
            m[i] = false;
            valid[i] = false;
        }
    }

    inline bool
    check(int x, int y) {
        y = y % wrow;
        switch (wrow) {
            case 3:
                return this->m[x] && 
                       this->m[x + width] && 
                       this->m[x + width * 2] &&
                       this->valid[x + y * width];     
            case 5:
                return this->m[x] && 
                       this->m[x + width] && 
                       this->m[x + width * 2] && 
                       this->m[x + width * 3] && 
                       this->m[x + width * 4] &&
                       this->valid[x + y * width];
            case 7:
                return this->m[x] && 
                       this->m[x + width] && 
                       this->m[x + width * 2] && 
                       this->m[x + width * 3] && 
                       this->m[x + width * 4] &&
                       this->m[x + width * 5] &&
                       this->m[x + width * 6] && 
                       this->valid[x + y * width];
            default:
                throw std::runtime_error("wrow must be one of 3, 5, 7.");               
        }
    }
};

template <typename T>
std::vector<int64_t>
processBlocksStratQueinnec(
    int numSamples,
    int numStrata,
    helper::RasterBandMetaData& band,
    access::Access& access,
    existing::Existing& existing,
    IndexStorageVectors& indices,
    IndexStorageVectors& queinnecIndices,
    FocalWindow& fw,
    std::vector<std::vector<OGRPoint>>& existingSamples,
    uint64_t multiplier,
    uint64_t queinnecMultiplier,
    xso::xoshiro_4x64_plus& rng,
    std::string allocation,
    OptimAllocationDataManager& optim,
    std::vector<double> weights,
    int width,
    int height) 
{
    T nanInt = static_cast<T>(band.nan);
 
    //adjust blocks to be a large chunk of scanlines
    int xBlockSize = band.xBlockSize;
    int yBlockSize = band.yBlockSize;
    if (xBlockSize != width) {
        xBlockSize = width;
    }
    else {
        yBlockSize = std::min(128, height);
    }
 
    //allocate required memory
    band.p_buffer = VSIMalloc3(xBlockSize, yBlockSize + fw.vpad * 2, band.size);
    T *p_buffer = reinterpret_cast<T *>(band.p_buffer);
    int8_t *p_access = nullptr;
    if (access.used) {
        access.band.p_buffer = VSIMalloc3(xBlockSize, yBlockSize, access.band.size);
        p_access = reinterpret_cast<int8_t *>(access.band.p_buffer);
    }
 
    helper::RandValController rand(xBlockSize, yBlockSize, multiplier, &rng);
    helper::RandValController queinnecRand(xBlockSize, yBlockSize, queinnecMultiplier, &rng);
 
    if (optim.used) {
        optim.init(numStrata, xBlockSize, yBlockSize);
    }
 
    //get number of blocks
    int yBlocks = (height + yBlockSize - 1) / yBlockSize;
 
    int fwyi = 0;
    for (int yBlock = 0; yBlock < yBlocks; yBlock++) {
        //read block
        int xOff = 0;
        int yOff = yBlock * yBlockSize;
        int xValid = width;
        int yValid = std::min(yBlockSize, height - yBlock * yBlockSize);
 
        //read block
        if (yBlock == 0) { 
            CPLErr err = band.p_band->RasterIO(GF_Read, xOff, yOff, xValid, yValid, band.p_buffer, xValid, yValid, band.type, 0, 0);
            if (err) {
                throw std::runtime_error("error reading block from raster.");
            }
        }
        else {
            int stratYOff = yOff - fw.vpad * 2;
            int stratYValid = yValid + fw.vpad * 2;
            CPLErr err = band.p_band->RasterIO(GF_Read, xOff, stratYOff, xValid, stratYValid, band.p_buffer, xValid, stratYValid, band.type, 0, 0);
            if (err) {
                throw std::runtime_error("error reading block from raster.");
            }
        }
 
        //read access block
        if (access.used) {
            CPLErr err = access.band.p_band->RasterIO(GF_Read, xOff, yOff, xValid, yValid, access.band.p_buffer, xValid, yValid, GDT_Int8, 0, 0);
            if (err) {
                throw std::runtime_error("error reading block from raster.");
            }
        }
 
        //read mraster block for optim
        if (optim.used) {
            CPLErr err = optim.band.p_band->RasterIO(GF_Read, xOff, yOff, xValid, yValid, optim.band.p_buffer, xValid, yValid, optim.band.type, 0, 0);
            if (err) {
                throw std::runtime_error("error reading block from raster.");
            }
        }
 
        //calculate rand vals
        rand.calculateRandValues();
        queinnecRand.calculateRandValues();
 
        //calculate the within-block index. In the case where there is a padding at the top of the block,
        //adjust for that.
        int newBlockStart = (yBlock == 0) ? 0 : width * fw.vpad * 2;
 
        //iterate through block and update vectors
        for (int y = 0; y < yValid; y++) {
            //reset upcomming row in focal window matrix
            fw.reset(yBlock * yBlockSize + y); 
            
            //the indexes from 0 to the horizontal pad cannot be a queinnec index
            //because they are too close to the edges of the raster 
            for (int x = 0; x < fw.hpad; x++) {
                T val = p_buffer[newBlockStart + y * width + x];
                helper::Index index = {x, y + yBlock * yBlockSize};
 
                bool isNan = val == nanInt;
                bool accessible = !access.used || p_access[y * width + x] != 1;
                bool alreadySampled = existing.used && existing.containsIndex(index.x, index.y);
 
                //check nan
                if (isNan) {
                    continue;
                }
 
                //update optim allocation variance calculations
                if (optim.used) {
                    optim.update(y * width + x, val);
                }
 
                //update strata counts
                indices.updateStrataCounts(val);
 
                //update existing samled strata
                if (alreadySampled) {
                    existingSamples[val].push_back(existing.getPoint(index.x, index.y));
                }
 
                //add val to stored indices
                if (accessible && !alreadySampled) {
                    indices.updateFirstXIndexesVector(val, index);
 
                    if (rand.next()) {
                        indices.updateIndexesVector(val, index);
                    }
                }
            }
    
            for (int x = fw.hpad; x < width - fw.hpad; x++) {
                T val = p_buffer[newBlockStart + y * width + x];
                helper::Index index = {x, y + yBlock * yBlockSize};
 
                bool isNan = val == nanInt;
                bool accessible = !access.used || p_access[y * width + x] != 1;
                bool alreadySampled = existing.used && existing.containsIndex(index.x, index.y);
                
                //check nan
                if (isNan) {
                    continue;
                }
 
                //update optim allocation variance calculations
                if (optim.used) {
                    optim.update(y * width + x, val);
                }
 
                //update strata counts
                indices.updateStrataCounts(val);
 
                //update exisitng sampled strata
                if (alreadySampled) {
                    existingSamples[val].push_back(existing.getPoint(index.x, index.y));
                };
 
                //set focal window matrix value by checking horizontal indices within focal window
                int start = newBlockStart + y * width + x - fw.hpad;
                switch (fw.wcol) {
                    case 3:
                        fw.m[fwyi + x] = p_buffer[start] == p_buffer[start + 1] &&
                                 p_buffer[start] == p_buffer[start + 2];
                        break;
                    case 5: 
                        fw.m[fwyi + x] = p_buffer[start] == p_buffer[start + 1] &&
                                 p_buffer[start] == p_buffer[start + 2] &&
                                 p_buffer[start] == p_buffer[start + 3] &&
                                 p_buffer[start] == p_buffer[start + 4];
                        break;
                    case 7: 
                        fw.m[fwyi + x] = p_buffer[start] == p_buffer[start + 1] &&
                                 p_buffer[start] == p_buffer[start + 2] &&
                                 p_buffer[start] == p_buffer[start + 3] &&
                                 p_buffer[start] == p_buffer[start + 4] &&
                                 p_buffer[start] == p_buffer[start + 5] &&
                                 p_buffer[start] == p_buffer[start + 6];
                        break;
                    default:
                        throw std::runtime_error("wcol must be one of 3, 5, 6.");
                }
                
                //add val to stored indices and update focal window validity    
                if (accessible && !alreadySampled) {
                    indices.updateFirstXIndexesVector(val, index);
 
                    if (rand.next()) {
                        indices.updateIndexesVector(val, index);
                    }
                    
                    fw.valid[fwyi + x] = true;
                }
 
                //add index to queinnec indices if surrounding focal window is all the same
                helper::Index fwIndex = {x, index.y - fw.vpad};
                if (fw.check(fwIndex.x, fwIndex.y)) {
                    bool add = false;
                    int start = newBlockStart + y * width + x - 2 * width * fw.vpad;
                    switch (fw.wrow) {
                        case 3:
                            add = p_buffer[start] == p_buffer[start + width * 1] &&
                                  p_buffer[start] == p_buffer[start + width * 2];
                            break;
                        case 5:
                            add = p_buffer[start] == p_buffer[start + width * 1] &&
                                  p_buffer[start] == p_buffer[start + width * 2] &&
                                  p_buffer[start] == p_buffer[start + width * 3] &&
                                  p_buffer[start] == p_buffer[start + width * 4];
                            break;
                        case 7:
                            add = p_buffer[start] == p_buffer[start + width * 1] &&
                                  p_buffer[start] == p_buffer[start + width * 2] &&
                                  p_buffer[start] == p_buffer[start + width * 3] &&
                                  p_buffer[start] == p_buffer[start + width * 4] &&
                                  p_buffer[start] == p_buffer[start + width * 5] &&
                                  p_buffer[start] == p_buffer[start + width * 6];
                            break;
                    }
 
                    if (add) {
                        //(we know if we've made it here that val is the same for both the fw add and the current index)
                        queinnecIndices.updateFirstXIndexesVector(val, fwIndex);
 
                        if (queinnecRand.next()) {
                            queinnecIndices.updateIndexesVector(val, fwIndex);
                        }
                    }
                }
            }
 
            //the indexes from width - horizontal pad to width cannot be a queinnec index
            //because they are too close to the edges of the raster 
            for (int x = width - fw.hpad; x < width; x++) {
                T val = p_buffer[newBlockStart + y * width + x];
                helper::Index index = {x, y + yBlock * yBlockSize};
 
                bool isNan = val == nanInt;
                bool accessible = !access.used || p_access[y * width + x] != 1;
                bool alreadySampled = existing.used && existing.containsIndex(index.x, index.y);
 
                //check nan
                if (isNan) {
                    continue;
                }
 
                //update optim allocation variance calculations
                if (optim.used) {
                    optim.update(y * width + x, val);
                }
 
                //update strata counts
                indices.updateStrataCounts(val);
 
                //update existing sampled strata
                if (alreadySampled) {
                    existingSamples[val].push_back(existing.getPoint(index.x, index.y));
                }
 
                //add val to stored indices
                if (accessible && !alreadySampled) {
                    indices.updateFirstXIndexesVector(val, index);
 
                    if (rand.next()) {
                        indices.updateIndexesVector(val, index);
                    }
                }
            }
 
            fwyi += width;
            if (fwyi == width * fw.wrow) {
                fwyi = 0;
            }
        }
    }
 
    VSIFree(band.p_buffer);
    if (access.used) {
        VSIFree(access.band.p_buffer);
    }
 
    if (optim.used) {
        weights = optim.getAllocationPercentages();
    }
    std::vector<int64_t> strataSampleCounts = calculateAllocation(
        numSamples,
        allocation,
        indices.getStrataCounts(),
        weights,
        indices.getNumDataPixels()
    );
 
    return strataSampleCounts;  
}

std::tuple<std::vector<std::vector<double>>, vector::GDALVectorWrapper *, size_t>
strat(
    raster::GDALRasterWrapper *p_raster,
    int bandNum,
    int64_t numSamples,
    int64_t numStrata,
    std::string allocation,
    std::vector<double> weights,
    raster::GDALRasterWrapper *p_mraster,
    int mrastBandNum,
    std::string method,
    int wrow,
    int wcol,
    double mindist,
    vector::GDALVectorWrapper *p_existing,
    bool force,
    vector::GDALVectorWrapper *p_access,
    std::string layerName,
    double buffInner,
    double buffOuter,
    std::vector<std::pair<std::string, int>> mapStratMapping,
    bool plot,
    std::optional<uint64_t> random_state,
    std::string filename,
    std::string tempFolder)
{
    GDALAllRegister();
 
    double mindist_sq = mindist * mindist;
    bool useMindist = mindist != 0;
    int width = p_raster->getWidth();
    int height = p_raster->getHeight();
    double *GT = p_raster->getGeotransform();
    bool mapped = mapStratMapping.size() != 0;      
 
    std::mutex bandMutex;
    std::mutex rngMutex;
    std::mutex accessMutex;
 
    //step 1: get raster band
    helper::RasterBandMetaData band;
 
    band.p_band = p_raster->getRasterBand(bandNum);
    band.type = p_raster->getRasterBandType(bandNum);
    band.size = p_raster->getRasterBandTypeSize(bandNum);
    band.p_buffer = nullptr;
    band.nan = band.p_band->GetNoDataValue();
    band.p_mutex = &bandMutex;
    band.p_band->GetBlockSize(&band.xBlockSize, &band.yBlockSize);
 
    helper::printTypeWarningsForInt32Conversion(band.type);
    
    //create output dataset before doing anything which will take a long time in case of failure.
    GDALDriver *p_driver = GetGDALDriverManager()->GetDriverByName("MEM");
    if (!p_driver) {
        throw std::runtime_error("unable to create output sample dataset driver.");
    }
    GDALDataset *p_samples = p_driver->Create("", 0, 0, 0, GDT_Unknown, nullptr);
    if (!p_samples) {
        throw std::runtime_error("unable to create output dataset with driver.");
    }
    
    vector::GDALVectorWrapper *p_wrapper = new vector::GDALVectorWrapper(p_samples, std::string(p_raster->getDataset()->GetProjectionRef()));
    OGRLayer *p_layer = p_samples->CreateLayer("samples", p_wrapper->getSRS(), wkbPoint, nullptr);
    if (!p_layer) {
        throw std::runtime_error("unable to create output dataset layer.");
    }
 
    if (!mapped) {
        OGRFieldDefn strataField(band.p_band->GetDescription(), OFTInteger);
        OGRErr err = p_layer->CreateField(&strataField);
        if (err) {
            throw std::runtime_error("cannot create field for strata name.");
        }
    }
    else {
        for (size_t i = 0; i < mapStratMapping.size(); i++) {
            OGRFieldDefn strataField(mapStratMapping[i].first.c_str(), OFTInteger);
            OGRErr err = p_layer->CreateField(&strataField);
            if (err) {
                throw std::runtime_error("cannot create create field for strata name.");
            }
        }
    }
 
    access::Access access(
        p_access, 
        p_raster, 
        layerName, 
        buffInner, 
        buffOuter, 
        true, 
        tempFolder, 
        band.xBlockSize,
        band.yBlockSize
    );
 
    std::vector<double> xCoords, yCoords;
    std::vector<std::vector<OGRPoint>> existingSamples(numStrata);  
    existing::Existing existing(
        p_existing,
        p_raster,
        GT,
        width,
        p_layer,
        false,
        xCoords,
        yCoords,
        false   
    );
 
    //fast random number generator using xoshiro256++
    //https://vigna.di.unimi.it/ftp/papers/ScrambledLinear.pdf
    xso::xoshiro_4x64_plus rng = random_state ? xso::xoshiro_4x64_plus(*random_state) : xso::xoshiro_4x64_plus();
    
    uint64_t multiplier = helper::getProbabilityMultiplier(
        width, 
        height, 
        p_raster->getPixelWidth(),
            p_raster->getPixelHeight(), 
        8, 
        numSamples, 
        useMindist, 
        access.area
    );
    uint64_t queinnecMultiplier = helper::getProbabilityMultiplier(
        width, 
        height, 
        p_raster->getPixelWidth(), 
        p_raster->getPixelHeight(), 
        64, 
        numSamples, 
        useMindist, 
        access.area
    );
 
    //normal indices used with both methods, queinnec indices used only with queinnec method
    IndexStorageVectors indices(numStrata, 10000);
    IndexStorageVectors queinnecIndices(numStrata, 10000);   
 
    FocalWindow fw(wrow, wcol, width);
 
    OptimAllocationDataManager optim(p_mraster, mrastBandNum, allocation);
 
    std::vector<int64_t> strataSampleCounts; 
    if (method == "random") {
        switch (band.type) {
            case GDT_Int8:
                strataSampleCounts = processBlocksStratRandom<int8_t>(numSamples, numStrata, band, access, 
                                              existing, indices, existingSamples,
                                              multiplier, rng, allocation, optim,
                                              weights, width, height);
                break;
            case GDT_Int16:
                strataSampleCounts = processBlocksStratRandom<int16_t>(numSamples, numStrata, band, access, 
                                              existing, indices, existingSamples,
                                              multiplier, rng, allocation, optim,
                                              weights, width, height);
                break;
            default:
                strataSampleCounts = processBlocksStratRandom<int32_t>(numSamples, numStrata, band, access, 
                                              existing, indices, existingSamples,
                                              multiplier, rng, allocation, optim,
                                              weights, width, height);
                break;
        }
    }
    else { //method == queinnec
        switch (band.type) {
            case GDT_Int8:
                strataSampleCounts = processBlocksStratQueinnec<int8_t>(numSamples, numStrata, band, access, existing, 
                                            indices, queinnecIndices, fw, existingSamples,
                                            multiplier, queinnecMultiplier, rng, allocation, 
                                            optim, weights, width, height);
                break;
            case GDT_Int16:
                strataSampleCounts = processBlocksStratQueinnec<int16_t>(numSamples, numStrata, band, access, existing, 
                                            indices, queinnecIndices, fw, existingSamples,
                                            multiplier, queinnecMultiplier, rng, allocation, 
                                            optim, weights, width, height);
                break;
            default:
                strataSampleCounts = processBlocksStratQueinnec<int32_t>(numSamples, numStrata, band, access, existing, 
                                            indices, queinnecIndices, fw, existingSamples,
                                            multiplier, queinnecMultiplier, rng, allocation, 
                                            optim, weights, width, height);
                break;
        }
    }
 
    std::vector<std::vector<helper::Index> *> strataIndexVectors;
    std::vector<size_t> nextIndexes(numStrata, 0);
 
    std::vector<bool> completedStrata(numStrata, false);
    std::vector<bool> completedStrataQueinnec(numStrata, false);
    std::vector<int64_t> samplesAddedPerStrata(numStrata, 0);
    int64_t numCompletedStrata = 0;
    int64_t numCompletedStrataQueinnec = 0;
    int64_t curStrata = 0;
    int64_t addedSamples = 0;
 
    //the following looks complicated. The reason for all of the fields, vectors, and pointers 
    //is to lower the ammount of times Fields and vectors are copied or created. Fields are
    //used to contain metadata on each sample, for example whether that sample is part of
    //an existing sample plot network, or the strata from which that field originated.
    //
    //PROBLEMS:
    //One option, would be to create a new Field or a vector of Fields EVERY iteration when a point
    //is added. Obviously, this incurs a large overhead.
    //
    //Another option would be to do without the pointers and have a set vector of fields or a set
    //vector of field vectors. However, whenever a vector is indexed it returns a copy of the element.
    //So, rather than creating a new Field or vector of fields, it would copy an existing field or
    //vector of fields EVERY iteration when a point is added. Obviously, this also incurs a large
    //overhead. 
    //
    //SOLUTION:
    //in order to get around these problems, ONLY the exact fields which will be used are ever 
    //created, which includes the 'fieldExistingTrue', 'fieldExistingFalse', and 'strataFields'
    //variables.
    //
    //THEN, since we may need to add multiple fields for each point, in addition to those fields
    //we also create every combination of fields that would ever be required in the form of
    //std::vector<Field *> (so the field is not copied). This is the 'fieldVectors' vector which
    //has the type std::vector<std::vector<Field *>>.
    //
    //Essentially: 'fieldExistingTrue', 'fieldExistingFalse', and 'strataFields' store the ONLY
    //copies of the fields (and are accessed via pointers.) 'fieldVectors' stores all the COMBINATIONS
    //of (pointers to) fields that might be relevant to a pixel (and are accessed via pointers.)
    //
    //Finally, to make the COMBINATIONS of fields easily indexible without copying, there are
    //additional vectors storing POINTERS to combinations within 'fieldVectors'. These are
    //'fieldVectorPointersExistingTrue', 'fieldVectorPointersExistingFalse', and 'fieldVectorPointersNoExisting',
    //and have the type std::vector<std::vector<Field *> *> in order to point to a particular
    //entry in the 'fieldVectors' vector.   
    helper::Field fieldExistingTrue("existing", 1);
    helper::Field fieldExistingFalse("existing", 0);
    std::vector<helper::Field> strataFields(mapped ? 0 : numStrata);
    std::vector<std::vector<helper::Field>> mappedStrataFields(mapped ? mapStratMapping.size() : 0);
    if (!mapped) {
        strataFields.resize(numStrata);
        for (int i = 0; i < numStrata; i++) {
            strataFields[i].fname = std::string(band.p_band->GetDescription());
            strataFields[i].fval = i;
        }
    }
    else { //mapped
        for (size_t i = 0; i < mapStratMapping.size(); i++) {
            mappedStrataFields[i].resize(mapStratMapping[i].second);
            for (int strat = 0; strat < mapStratMapping[i].second; strat++) {
                mappedStrataFields[i][strat].fname = mapStratMapping[i].first;
                mappedStrataFields[i][strat].fval = strat;
            }
        }
    }
    
    int fieldVectorCount = existing.used ? numStrata * 2 : numStrata;
    std::vector<std::vector<helper::Field *>> fieldVectors(fieldVectorCount);
    
    int fvi = 0; //field vectors index
    std::vector<std::vector<helper::Field *> *> fieldVectorPointersExistingTrue(existing.used ? numStrata : 0);
    std::vector<std::vector<helper::Field *> *> fieldVectorPointersExistingFalse(existing.used ? numStrata : 0);
    std::vector<std::vector<helper::Field *> *> fieldVectorPointersNoExisting(existing.used ? 0 : numStrata);
 
    if (existing.used && !mapped) {
        for (int strata = 0; strata < numStrata; strata++) {
            //resize entries for both existing==true and existing==false
            fieldVectors[fvi].resize(2);
            fieldVectors[fvi + 1].resize(2);
 
            //update field vector pointers for strata and existing
            fieldVectors[fvi][0] = &fieldExistingTrue;
            fieldVectors[fvi + 1][0] = &fieldExistingFalse;
            fieldVectors[fvi][1] = strataFields.data() + strata;
            fieldVectors[fvi + 1][1] = strataFields.data() + strata;
 
            //update pointers to field vectors
            fieldVectorPointersExistingTrue[strata] = fieldVectors.data() + fvi;
            fieldVectorPointersExistingFalse[strata] = fieldVectors.data() + fvi + 1;
 
            //increment field vector index
            fvi += 2;
        }
    }
    else if (existing.used && mapped) {
        for (int i = 0; i < numStrata; i++) {
            int strata = i;
 
            //resize entries for both existing==true and existing==false
            fieldVectors[fvi].resize(mapStratMapping.size() + 1);
            fieldVectors[fvi + 1].resize(mapStratMapping.size() + 1);
 
            //update field vector pointers for strata and existing
            for (size_t j = 0; j < mapStratMapping.size(); j++) {
                int bandStrata = strata % mapStratMapping[j].second;
                fieldVectors[fvi][j] = mappedStrataFields[j].data() + bandStrata;
                fieldVectors[fvi + 1][j] = mappedStrataFields[j].data() + bandStrata;
                strata = strata / mapStratMapping[j].second;
            }
            fieldVectors[fvi][mapStratMapping.size()] = &fieldExistingTrue;
            fieldVectors[fvi + 1][mapStratMapping.size()] = &fieldExistingFalse;
 
            //update pointers to field vectors
            fieldVectorPointersExistingTrue[i] = fieldVectors.data() + fvi;
            fieldVectorPointersExistingFalse[i] = fieldVectors.data() + fvi + 1;
 
            //increment field vector index
            fvi += 2;
        }   
    }
    else if (!existing.used && !mapped) {
        for (int strata = 0; strata < numStrata; strata++) {
            fieldVectors[fvi].resize(1);
            fieldVectors[fvi][0] = strataFields.data() + strata;
            fieldVectorPointersNoExisting[strata] = fieldVectors.data() + fvi;
            
            fvi++;
        }
    }
    else { 
        for (int i = 0; i < numStrata; i++) {
            int strata = i;
            
            fieldVectors[fvi].resize(mapStratMapping.size());
            for (size_t j = 0; j < mapStratMapping.size(); j++) {
                int bandStrata = strata % mapStratMapping[j].second;
                fieldVectors[fvi][j] = mappedStrataFields[j].data() + bandStrata;
                strata = strata / mapStratMapping[j].second;
            }
            fieldVectorPointersNoExisting[i] = fieldVectors.data() + fvi;
            fvi++;
        }
    }
 
    //add existing sample plots
    helper::NeighborMap neighbor_map;
    if (existing.used) {
        if (force) {
            //if force is used, add all samples no matter what
            for (int i = 0; i < numStrata; i++) {
                std::vector<OGRPoint> samples = existingSamples[i];
                    for (const OGRPoint& point : samples) {
                    helper::addPoint(&point, p_layer, fieldVectorPointersExistingTrue[i]);
 
                    addedSamples++;
                    samplesAddedPerStrata[i]++;
 
                    if (plot) {
                        xCoords.push_back(point.getX());
                        yCoords.push_back(point.getY());
                    }
                }   
 
                if (samplesAddedPerStrata[i] >= strataSampleCounts[i]) {
                    numCompletedStrata++;
                    numCompletedStrataQueinnec++;
                    completedStrata[i] = true;
                    completedStrataQueinnec[i] = true;
                }
            }
        }
        else { //force == false
            for (size_t i = 0; i < numStrata; i++) {
                std::vector<OGRPoint> samples = existingSamples[i];
                std::shuffle(samples.begin(), samples.end(), rng);
                
                size_t j = 0;
                while (samplesAddedPerStrata[i] < strataSampleCounts[i] && j < samples.size()) {
                    OGRPoint point = samples[j];
                    double x = point.getX();
                    double y = point.getY();
                    bool valid = true;
 
                    if (useMindist) {
                        valid = helper::is_valid_sample(x, y, neighbor_map, mindist, mindist_sq);
                    }
 
                    if (valid) { 
                        helper::addPoint(&point, p_layer, fieldVectorPointersExistingTrue[i]);
 
                        addedSamples++;
                        samplesAddedPerStrata[i]++;
 
                        if (plot) {
                            xCoords.push_back(point.getX());
                            yCoords.push_back(point.getY());
                        }
                    }
 
                    j++;
                }
 
                if (samplesAddedPerStrata[i] == strataSampleCounts[i]) {
                    numCompletedStrata++;
                    numCompletedStrataQueinnec++;
                    completedStrata[i] = true;
                    completedStrataQueinnec[i] = true;
                }
            }
        }
    }   
 
    if (method == "Queinnec") {
        strataIndexVectors = queinnecIndices.getStrataIndexVectors(samplesAddedPerStrata, strataSampleCounts, rng);
 
        int64_t curStrata = 0;
        while (numCompletedStrataQueinnec < numStrata && addedSamples < numSamples) {
            if (curStrata == numStrata) {
                curStrata = 0;
            }
            if (completedStrataQueinnec[curStrata]) {
                curStrata++;
                continue;
            }
 
            int64_t sampleCount = strataSampleCounts[curStrata];
            int64_t samplesAdded = samplesAddedPerStrata[curStrata];
            if (samplesAdded == sampleCount) {
                numCompletedStrataQueinnec++;
                completedStrataQueinnec[curStrata] = true;
                numCompletedStrata++;
                completedStrata[curStrata] = true;
                curStrata++;
                continue;
            }
 
            std::vector<helper::Index> *strataIndexes = strataIndexVectors[curStrata];
            size_t nextIndex = nextIndexes[curStrata];  
            if (strataIndexes->size() == nextIndex) {
                numCompletedStrataQueinnec++;
                completedStrataQueinnec[curStrata] = true;
                curStrata++;
                continue;
            }
 
            helper::Index index = strataIndexes->at(nextIndex);
            nextIndexes[curStrata]++;
 
            bool valid = true;
            const auto [x, y] = helper::sample_to_point(GT, index);
            OGRPoint newPoint = OGRPoint(x, y);
 
            if (useMindist) {
                valid = helper::is_valid_sample(x, y, neighbor_map, mindist, mindist_sq);
            }
            
            if (valid) { 
                    existing.used ? 
                    helper::addPoint(&newPoint, p_layer, fieldVectorPointersExistingFalse[curStrata]) :
                    helper::addPoint(&newPoint, p_layer, fieldVectorPointersNoExisting[curStrata]);
                
                addedSamples++;
                samplesAddedPerStrata[curStrata]++;
                
                if (plot) {
                    xCoords.push_back(x);
                    yCoords.push_back(y);
                }
            }
 
            curStrata++;
        }
    }
        
    strataIndexVectors = indices.getStrataIndexVectors(samplesAddedPerStrata, strataSampleCounts, rng);    
    for (int64_t i = 0; i < numStrata; i++) {
        //set next indexes to 0 because the vectors are different than the queinnec vectors
        nextIndexes[i] = 0;
    }
    curStrata = 0;
    
    //step 8: generate coordinate points for each sample index.
    while (numCompletedStrata < numStrata && addedSamples < numSamples) {
        if (curStrata == numStrata) {
            curStrata = 0;
        }
        if (completedStrata[curStrata]) {
            curStrata++;
            continue;
        }
 
        int64_t sampleCount = strataSampleCounts[curStrata];
        int64_t samplesAdded = samplesAddedPerStrata[curStrata];
        if (samplesAdded == sampleCount) {
            numCompletedStrata++;
            completedStrata[curStrata] = true;
            curStrata++;
            continue;
        }
 
        std::vector<helper::Index> *strataIndexes = strataIndexVectors[curStrata];
        size_t nextIndex = nextIndexes[curStrata];
        if (strataIndexes->size() == nextIndex) {
            numCompletedStrata++;
            completedStrata[curStrata] = true;
            curStrata++;
            continue;
        }
 
        helper::Index index = strataIndexes->at(nextIndex);
        nextIndexes[curStrata]++;
        
        bool valid = true;
        const auto [x, y] = helper::sample_to_point(GT, index);
        OGRPoint newPoint = OGRPoint(x, y);
 
        if (useMindist) {
            valid = helper::is_valid_sample(x, y, neighbor_map, mindist, mindist_sq);
        }
 
        if (valid) {
            existing.used ? 
                helper::addPoint(&newPoint, p_layer, fieldVectorPointersExistingFalse[curStrata]) :
                helper::addPoint(&newPoint, p_layer, fieldVectorPointersNoExisting[curStrata]);
                
        
            addedSamples++;
            samplesAddedPerStrata[curStrata]++;
 
            if (plot) {
                xCoords.push_back(x);
                yCoords.push_back(y);
            }
        }
 
        curStrata++;
    }
    
    //step 10: write vector if filename is not "".
    if (filename != "") {
        try {
            p_wrapper->write(filename);
        }
        catch (const std::exception& e) {
            std::cout << "Exception thrown trying to write file: " << e.what() << std::endl;
        }
    }
 
    size_t actualSampleCount = static_cast<size_t>(p_layer->GetFeatureCount());
    return {{xCoords, yCoords}, p_wrapper, actualSampleCount};
}
 
} //namespace strat
} //namespace sgs