Files @ 82182a2f068a
Branch filter:

Location: AENC/resampling_chain/montecarlo/resampler_cycle.cpp

Tom Bannink
Initial commit with Monte Carlo sampler
#include <iostream>
#include <random>
#include "resampler.hpp"

std::mt19937 mt(std::random_device{}());

class CycleResampler : public ResamplerBase<CycleResampler, std::mt19937> {
  public:
    CycleResampler(float p_, size_t n, bool allBads = true)
        : ResamplerBase(p_, n, mt) {
        assert(n >= 3);
        if (allBads) {
            for (auto i = 0u; i < n; ++i) {
                setVertex(i, BAD);
            }
        } else {
            setVertex(0, BAD);
        }
    }

    // Cycle has periodic boundary conditions
    void resampleNeighbors(int i) {
        resampleVertex(i);

        if (i == 0)
            resampleVertex(getN() - 1);
        else
            resampleVertex(i - 1);

        if (i + 1 == (int)getN())
            resampleVertex(0);
        else
            resampleVertex(i + 1);
    }

    // Returns number of steps before all good
    int run(int limit = -1) {
        int t;
        for (t = 0; t != limit && numBads() != 0; t++)
            doMove();
        return t;
    }
};

int main() {
    //if (0) {
        CycleResampler sampler(0.6f, 40);
        int loc;
        for (int i = 0; i < 2000; ++i) {
            loc = sampler.doMove();
            if (loc == -1)
                break;
            std::cout << '|';
            int j = 0;
            for (auto x : sampler.getState()) {
                if (j++ == loc) {
                    std::cout << 'X';
                } else {
                    std::cout << (x == BAD ? '*' : ' ');
                }
            }
            std::cout << '|' << std::endl;
        }
    //}

    // Measure mixing time for various n,p
    // starting from all-1 state
    if (0) {
        std::cout << "{{0,0,0}";
        for (int i = 0; i <= 8; ++i) {
            float p = 0.6f + 0.01f * i;
            for (int n = 5; n <= 25; ++n) {
                std::cerr << "Running p = " << p << " and n = " << n
                          << std::endl;
                for (int j = 0; j < 1000; ++j) {
                    CycleResampler sampler(p, n);
                    int time = sampler.run(500000);
                    std::cout << ',' << '{' << p << ',' << n << ',' << time
                              << '}';
                }
                std::cout << std::flush;
            }
        }
        std::cout << '}' << std::endl;
    }

    // Measure probability of hitting all-good
    // starting from single-bad and stopping if some condition is met
    // Interesting range is p in [0.60,0.75]
    std::cout << "{{0,0,0}";
    for (int i = 0; i <= 100; ++i) {
        float p = 0.60f + 0.001f * i;
        for (size_t n = 1000; n <= 1000; n += 200) {
            std::cerr << "Running p = " << p << " and n = " << n << std::endl;
            constexpr int samples = 500;
            int hitCount = 0;
            int otherHitCount = 0;
            for (int j = 0; j < samples; ++j) {
                CycleResampler sampler(p, n, false);
                for (int t = 0; t != 500000; t++) {
                    sampler.doMove();
                    if (sampler.numBads() == 0) {
                        hitCount++;
                        break;
                    }
                    //if (sampler.numBads() >= (5*n)/10) {
                    if (sampler.getState()[n/2] == BAD) {
                        otherHitCount++;
                        break;
                    }
                }
            }
            std::cout << ',' << '{' << p << ',' << n << ','
                      << float(hitCount) / float(samples) << '}';
            std::cout << std::flush;

            std::cerr << (samples - otherHitCount - hitCount) << '/' << samples
                      << " runs with timelimit.\n";
        }
    }
    std::cout << '}' << std::endl;

    return 0;
}