LatinHyperCube.h 3.97 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
//
// Class LatinHyperCube
//   This class does Latin hypercube sampling.
//
// Copyright (c) 2018, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// Implemented as part of the PhD thesis
// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
//
// This file is part of OPAL.
//
// OPAL 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.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
frey_m's avatar
frey_m committed
21 22 23 24 25 26
#ifndef OPAL_LATIN_HYPERCUBE_H
#define OPAL_LATIN_HYPERCUBE_H

#include "Sample/SamplingMethod.h"
#include "Sample/RNGStream.h"

frey_m's avatar
frey_m committed
27
#include <algorithm>
frey_m's avatar
frey_m committed
28
#include <deque>
frey_m's avatar
frey_m committed
29

frey_m's avatar
frey_m committed
30 31 32 33 34 35
class LatinHyperCube : public SamplingMethod
{

public:
    typedef typename std::uniform_real_distribution<double> dist_t;

frey_m's avatar
frey_m committed
36 37 38
    LatinHyperCube(double lower, double upper)
        : binsize_m(0.0)
        , upper_m(upper)
frey_m's avatar
frey_m committed
39
        , lower_m(lower)
frey_m's avatar
frey_m committed
40
        , dist_m(0.0, 1.0)
41 42
        , RNGInstance_m(RNGStream::getInstance())
        , seed_m(RNGStream::getGlobalSeed())
43
    {}
frey_m's avatar
frey_m committed
44
    
frey_m's avatar
frey_m committed
45 46 47
    LatinHyperCube(double lower, double upper, int seed)
        : binsize_m(0.0)
        , upper_m(upper)
frey_m's avatar
frey_m committed
48
        , lower_m(lower)
frey_m's avatar
frey_m committed
49
        , dist_m(0.0, 1.0)
50
        , RNGInstance_m(nullptr)
51
        , seed_m(seed)
52
    {}
frey_m's avatar
frey_m committed
53 54

    ~LatinHyperCube() {
55 56
        if ( RNGInstance_m )
            RNGStream::deleteInstance(RNGInstance_m);
frey_m's avatar
frey_m committed
57 58 59 60 61 62 63 64 65
    }

    void create(boost::shared_ptr<SampleIndividual>& ind, std::size_t i) {
        /* values are created within [0, 1], thus, they need to be mapped
         * the domain [lower, upper]
         */
        ind->genes[i] = map2domain_m(RNGInstance_m->getNext(dist_m));
    }
    
66 67 68 69 70 71 72 73 74 75 76
    void allocate(const CmdArguments_t& args, const Comm::Bundle_t& comm) {
        int id = comm.island_id;
        
        if ( !RNGInstance_m )
            RNGInstance_m = RNGStream::getInstance(seed_m + id);
        
        int nSamples = args->getArg<int>("nsamples", true);
        int nMasters = args->getArg<int>("num-masters", true);
        
        int nLocSamples = nSamples / nMasters;
        int rest = nSamples - nMasters * nLocSamples;
77
        
78 79
        if ( id < rest )
            nLocSamples++;
frey_m's avatar
frey_m committed
80
        
81
        int startBin = 0;
frey_m's avatar
frey_m committed
82
        
83 84 85 86 87 88 89 90 91 92 93 94
        if ( rest == 0 )
            startBin = nLocSamples * id;
        else {
            if ( id < rest ) {
                startBin = nLocSamples * id;
            } else {
                startBin = (nLocSamples + 1) * rest + (id - rest) * nLocSamples;
            }
        }
        
        binsize_m = ( upper_m - lower_m ) / double(nSamples);
        
95
        this->fillBins_m(nSamples, nLocSamples, startBin, seed_m);
frey_m's avatar
frey_m committed
96 97
    }
    
frey_m's avatar
frey_m committed
98 99 100 101 102 103 104 105 106 107 108 109
private:
    double map2domain_m(double val) {
        /* y = mx + q
         * 
         * [0, 1] --> [a, b]
         * 
         * y = (b - a) * x + a
         * 
         * where a and b are the lower, respectively, upper
         * bound of the current bin.
         */
        
frey_m's avatar
frey_m committed
110 111 112
        std::size_t bin = bin_m.back();
        bin_m.pop_back();
        
frey_m's avatar
frey_m committed
113
        return  binsize_m * (val + bin) + lower_m;
frey_m's avatar
frey_m committed
114
    }
frey_m's avatar
frey_m committed
115
    
116 117 118 119 120 121 122 123
    void fillBins_m(std::size_t nTotal, std::size_t nLocal, int startBin,
                    std::size_t seed)
    {
        std::deque<std::size_t> tmp;
        tmp.resize(nTotal);
        std::iota(tmp.begin(), tmp.end(), 0);

        // all masters need to shuffle the same way
124
        std::mt19937_64 eng(seed);
125 126 127 128 129 130
        std::shuffle(tmp.begin(), tmp.end(), eng);

        // each master takes its bins
        std::copy(tmp.begin()+startBin,
                  tmp.begin()+startBin+nLocal,
                  std::back_inserter(bin_m));
frey_m's avatar
frey_m committed
131
    }
frey_m's avatar
frey_m committed
132 133

private:
frey_m's avatar
frey_m committed
134
    std::deque<std::size_t> bin_m;
frey_m's avatar
frey_m committed
135 136
    double binsize_m;
    
frey_m's avatar
frey_m committed
137
    double upper_m;
frey_m's avatar
frey_m committed
138 139
    double lower_m;
    
frey_m's avatar
frey_m committed
140
    dist_t dist_m;
141 142
    
    RNGStream *RNGInstance_m;
143 144
    
    std::size_t seed_m;
frey_m's avatar
frey_m committed
145 146 147
};

#endif