sampling.h

Go to the documentation of this file.

/*
 *  helper/sampling.h - statistical sampling functions
 *  Copyright (C) 2015, D Haley 

 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.

 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.

 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef ATOMPROBE_SAMPLING_H
#define ATOMPROBE_SAMPLING_H
#include <vector>

#include "atomprobe/primitives/ionHit.h"
#include "atomprobe/atomprobe.h"
#include "maths/lfsr.h"

#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>

namespace AtomProbe 
{
    //TODO : Deprecate me
    //Randomly sample, without replacement. 
    void sampleIons(const std::vector<IonHit> &ions, float sampleFactor, 
                std::vector<IonHit> &sampled, bool strongRandom=true);

    //Randomly select subset. Subset ordering is not guaranteed to be either random or ordered 
    // Returns -1 on abort, otherwise returns number of randomly selected items
    // if the supplied RNG is null, a "fast" linear shift algorithm, with bad randomisation will be used
    template<class T> size_t randomSelect(std::vector<T> &result, 
                const std::vector<T> &source, size_t num,gsl_rng *rng)
    {
        //If there are not enough points, just copy it across in whole
        if(source.size() <= num)
        {
            num=source.size();
            result.resize(source.size());
#pragma omp parallel for
            for(size_t ui=0; ui<num; ui++)
                result[ui] = source[ui]; 
        
            return num;
        }

        result.resize(num);

        if(rng)
        {
            result.resize(num);

            gsl_ran_choose (rng, &(result[0]), num, (void*)&(source[0]), source.size(), sizeof(IonHit));

        }   
        else
        {
            //Use a weak randomisation
            LinearFeedbackShiftReg l;

            //work out the mask level we need to use
            size_t i=1;
            unsigned int j=0;
            while(i < (source.size()<<1))
            {
                i=i<<1;
                j++;
            }

            //linear shift table starts at 3.
            if(j<3) {
                j=3;
                i = 1 << j;
            }

            size_t start;
            //start at a random position  in the linear state
            start =(size_t)((double)rand()/RAND_MAX*i);
            l.setMaskPeriod(j);
            l.setState(start);

            size_t ui=0;    
            //generate unique weak random numbers.
            while(ui<num)
            {
                size_t res;
                res= l.clock();
                
                //use the source if it lies within range.
                //drop it silently if it is out of range
                if(res< source.size())
                {
                    result[ui] =source[res];
                    ui++;
                }
            }

        }

        return num;
    }

    //This function will not be efficient if num is close to nMax, and nMax is alrge
    // Should have inversion detection to improve efficiency
    template<class T> void randomIndices(std::vector<T> &res, size_t num, size_t nMax, gsl_rng *rng)
    {
        num = std::min(num,nMax);
        res.resize(num);

        //Use reservoir sampling
        for(auto i=0; i<num; i++)
            res[i]=i;

        for(auto i=num; i<nMax;i++)
        {
            unsigned int j;
            j = gsl_rng_uniform_int(rng,i);

            if(j <num)
                res[j]=i;

        }
        
    }
    


#ifdef DEBUG
//Unit tests
bool testSampling();
#endif
}
#endif