diff --git a/montecarlo/resampler_segment.cpp b/montecarlo/resampler_segment.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6ffd1b74b21564cecdb1ceafa091c03566092e67 --- /dev/null +++ b/montecarlo/resampler_segment.cpp @@ -0,0 +1,92 @@ +#include +#include +#include "resampler.hpp" + +std::mt19937 mt(std::random_device{}()); + +class SegmentResampler : public ResamplerBase { + public: + SegmentResampler(float p_, size_t n) : ResamplerBase(p_, n, mt) { + setVertex(0, BAD); + } + + // No periodic boundary + void resampleNeighbors(int i) { + resampleVertex(i); + + if (i != 0) + resampleVertex(i - 1); + if (i + 1 != (int)getN()) + resampleVertex(i + 1); + } +}; + +void doRun(float p, size_t n) { + std::cerr << "Running p = " << p << " and n = " << n << std::endl; + + constexpr int samples = 1000; + + int terminateCount = 0; + int hitnCount = 0; + for (int j = 0; j < samples; ++j) { + SegmentResampler sampler(p, n); + for (auto t = 0u; t != 1000 * n; t++) { + sampler.doMove(); + if (sampler.numBads() == 0) { + terminateCount++; + break; + } + if (sampler.getState()[n - 1] == BAD) { + hitnCount++; + break; + } + } + } + if (terminateCount + hitnCount != samples) { + std::cerr << (samples - terminateCount - hitnCount) << '/' << samples + << " runs with timelimit!\n"; + } + + std::cout << '{' << p << ',' << n; + std::cout << ',' << hitnCount << '/' << hitnCount + terminateCount; + std::cout << ',' << hitnCount + terminateCount; + std::cout << ',' << samples; + std::cout << '}'; + std::cout << std::flush; +} + +int main() { + // Measure probability of hitting vertex N starting from single-bad + // Interesting range is p in [0.60,0.75] + std::cout << "(* Probability of hitting site n when starting from single " + "BAD at 0 on segment [0,n].\n"; + std::cout << "Data-format: {p, n, probability, runs within timelimit, " + "total run attempts} *)\n"; + std::cout << "{\n"; + for (size_t n = 600;; n += 200) { + doRun(0.10f, n); + std::cout << ',' << std::endl; + + for (int i = 0; i <= 100; ++i) { + float p = 0.60f + 0.001f * i; + doRun(p, n); + std::cout << ',' << std::endl; + } + doRun(0.75f, n); + std::cout << ',' << std::endl; + ; + doRun(0.80f, n); + std::cout << ',' << std::endl; + ; + doRun(0.90f, n); + + if (n == 1000) + break; + else + std::cout << ',' << std::endl; + ; + } + std::cout << '\n' << '}' << std::endl; + + return 0; +}