libcrbn/sampler/latinHypercubeSampler.hpp

#ifndef _LATIN_HYPERCUBE_SAMPLER_HPP_
#define _LATIN_HYPERCUBE_SAMPLER_HPP_

#include <string.h>
#include "sampler.h"

// TODO : add home random function

// Latin hypercube sampling
class LatinHypercubeSampler : public Sampler {
  public:
        // Constructor
        LatinHypercubeSampler(const int& iNSamples,
                              const int& iNDimensions) :
                Sampler(iNSamples, iNDimensions)
        {
                _samples = new float[_nSamples * _nDimensions];
                _samplePos = _nSamples;
        }

        // Destructor
        virtual
        ~LatinHypercubeSampler() {
                if(_samples) delete [] _samples;
                _samples = 0;
        }
        
        // Computes next sample and store it in the array
        // passed as argument.
        // Returns false if we generated the s'th sample.
        bool next(float* iV) {
                if(_samplePos == _nSamples) {
                        init();
                        _samplePos = 0;
                }
                
                memcpy (iV, _samples + (_samplePos * _nDimensions), sizeof(float[_nDimensions]));
                ++_samplePos;
                return (_samplesPos < _nSamples);
        }

  private:
        // Initialize samples array
        void init() {
                int i, j;
                float delta;

                // Generate LHS samples along the diagonal
                delta = 1.0 / _nSamples;

                for(j = 0; j < _nSamples; ++j) {
                        for(i = 0; i < _nDimensions; ++i) {
                                _samples[i + (j * _nDimensions)] = (j + (random() / (float)RAND_MAX)) * delta;
                        }
                }

                // Permute LHS in each dimension
                for(i = 0; i < _nDimensions; ++i) {
                        for(j = 0; j < _nSamples; ++j) {
                                int other = random() % _nSamples;
                                float tmp = _samples[i + (j * _nDimensions)];

                                _samples[i + (j     * nDimensions)] =  _samples[i + (other * _nDimensions)];
                                _samples[i + (other * _nDimensions)] =  tmp;
                        }
                }
        }
        
  private:
        // Samples
        float* _samples;
};

#endif

---