ivajankovic/shake_and_bake_report.md

## shake_and_bake_report.md

      
    Raw
  

              shake_and_bake_report.md
            
          
    GSOC Final Submission Report

Project: Shake and Bake - Sampling from the boundary of convex polytopes
Contributor: Iva Janković
Mentors: Vissarion Fisikopoulos, Elias Tsigaridas, Apostolos Chalkis (GeomScale)
Overview

The main focus of the project was to implement the Markov Chain Monte Carlo-based Shake and Bake (SB) algorithms for uniform sampling from the boundary of a convex polytope and integrate it into the C++ library volesti, along with its R interface Rvolesti and the Python package dingo. In addition, another deliverable was added to the project, a new variant, the Billiard Shake-and-Bake, that has been implemented and tested. Finally, what is left to do, is test the functionality of the implemented algorithms on several metabolic networks of varying dimensionalities, where convergence performance will be evaluated using PSRF and ESS diagnostics.
Shake and Bake Method

The SB method generates a Markov Chain whose stationary distribution converges to the uniform distribution over the boundary of a convex polytope defined by a set of linear inequality constraints. It performs a random walk on ∂S₀ as follows:

Initialization – choose a starting boundary point.  - compute_boundary_point()
Direction Sampling – generate a random direction vector.  - get_direction() within initialize()
Compute Hit Point – find the intersection of the direction with the polytope boundary.  - apply()
Move Decision – decide whether to accept the new point (Running SB: acceptance = 1 ). - apply()
Repeat – continue until the stopping criterion is met.

It follows the steps as described in [1]. To test the uniformity, the scaling ratio test along with random facet 2D projection for generic polytopes (cube, simplex, birkhoff) with the plots are added in sbtest.py (with generated txt files).
Original /Limping Shake and Bake walk

The implemented SB variant is Running SB because it performs significantly better than its counterparts Original and Limping. But, if one wants to test that out, here is the additional piece of code to be added in apply function alongside some simple enum switch / branch logic. Also, it is very important to use  _v = GetDirection<Point>::apply(n, rng);  (coressponds to Step 1. of Original and Limping SB from [1] ) instead of Point v = get_direction(P,rng); (for Running SB).
NT beta;
if (mode_ == Mode::Original) {
    NT den = dot_r - dot_k;
    if (std::abs(den) < eps) continue;
    beta = std::clamp(dot_r / den, NT(0), NT(1));
} 
else {
    beta = -dot_k;
}

if (beta > NT(0) && beta <= NT(1) &&
    rng.sample_urdist() < beta)
{
    p_         = y;
    facet_idx_ = facet_new;
    A_row_k_   = A_row_r;
    Ar_.noalias() -= lambda_hit * Av_;  
}
Billiard Shake and Bake walk

This is a variant of Running SB algorithm implemented in that also belongs to the class of uniform boundary sampling algorithms. It follows the steps as described in [1], but after Step 1 (direction sampling) it does number of reflections ( with compute_reflection) defined by upper bound of reflection (nr). For each step, the number of reflections is sampled by using inverse exponential distribution, i.e. the number of reflections is (1-z)*nr. There is also an option for sampling uniformly from {1,...,nr} that is commented out.
Volesti / Rvolesti / dingo updates

PR for volesti SB implementation GeomScale/volesti#361
PR for volesti Billiard SB implementation GeomScale/volesti#365
The layout is as follows:

The algorithm logic - random_walks/shake_and_bake_walk.hpp / random_walks/billiard_shake_and_bake_walk.hpp
An example code of usage and testing - examples/sb_sampling / examples/bsb_sampling
The unit tests - tests/shake_and_bake_test.cpp / tests/billiard_shake_and_bake_test.cpp

Changes and that were made in order for the implementation to work:


The implementation allows for the user to give the starting boundary point, but in case it is not provided, a function for finding a boundary point of a given polytope called compute_boundary_point(P,rng, eps) was created. (File: preprocess/feasible_point.hpp)


Since these algorithms heavily rely on the facet indexing, the existing functions (line_positive_intersect and similar) for finding the intersection were modify to accept additional parameter called facet_idx as a facet that will be skipped in the search. (File: convex_bodies/hpolytope.h)


To shorten the code a function get_row(facet_idx) that returns a normal vector of the facet that a boundary point lies on was implemented. (File: convex_bodies/hpolytope.h)


The scaling_ratio_boundary_test function takes a set of boundary samples from a polytope and, for each facet, finds the samples lying on it, computes their centroid, and then repeatedly rescales the polytope around that centroid. For each rescaling step, it measures how many of the original facet samples remain feasible in the shrunken polytope, and compares these survival ratios with the theoretical scaling law x^(dim). By recording the deviations (average and maximum) across scaling factors, the test shows how uniformly the sampling covers each facet and whether some facets are underrepresented in the boundary distribution. (File: diagnostics/scaling_ratio.hpp)


In Billiard SB, for sampling a number of reflections per step truncated exponential is used with help of a new function sample_trunc_expdist() that generates a sample from an exponential distribution until the value is ≤ 1. (File: generators/boost_random_number_generator.hpp)


Structs BilliardShakeAndBakeRandomPointGenerator and ShakeAndBakeRandomPointGenerator were implemented with coressponding sampling functions billiard_shakeandbake_sampling and shakeandbake_sampling for unit tests, rvolesti and dingo integration. (Files: sampling/sampling.hpp & random_point_generators.hpp)


PR for rvolesti implementations GeomScale/Rvolesti#34
Rcpp wrappers for both SB and Billiard SB are in sample_points.cpp. Example of usage:
P = gen_cube(3, 'H')
points = sample_points(P, n = 100, random_walk = list("walk" = "SB", "walk_length" = 5))
points = sample_points(P, n = 100, random_walk = list("walk" = "BSB", "walk_length" = 5))
In the continuation of this work, our Rcpp wrappers will be compared with existing R implementation od SB [2], alongside other boundary sampling algorithms from volesti.
PR for dingo implementations  GeomScale/dingo#119
References

[1] C. G. E. Boender, R. J. Caron, J. F. McDonald, A. H. G. Rinnooy Kan, H. E. Romeijn, R. L. Smith, J. Telgen i A. C. F. Vorst,

Shake-And-Bake Algorithms for Generating Uniform Points on the Boundary of Bounded Polyhedra, 1991.

Available at: https://doi.org/10.1287/opre.39.6.945
[2] "Shake and Bake" sampler, https://rdrr.io/cran/hitandrun/man/shakeandbake.html.
No results found