Implementation of the Shake and Bake algorithm and the Scaling Ratio diagnostic

include/convex_bodies/hpolytope.h

@@ -7,6 +7,7 @@
 //Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program.
 //Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
 //Contributed and/or modified by Luca Perju, as part of Google Summer of Code 2024 program.
+//Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program.
 
 // Licensed under GNU LGPL.3, see LICENCE file
 
@@ -203,6 +204,13 @@ class HPolytope {
         return b;
     }
 
+    VT get_row(int facet_index) const 
+    {
+        if (facet_index < 0 || facet_index >= A.rows())
+            throw std::out_of_range("Facet index out of range!");
+        return A.row(facet_index);
+    }
+
     bool is_normalized ()
     {
         return normalized;
@@ -535,6 +543,7 @@ class HPolytope {
         NT inner_prev = params.inner_vi_ak;
         VT sum_nom;
         int m = num_of_hyperplanes(), facet;
+        int skip = params.facet_prev;
 
         Ar.noalias() += lambda_prev*Av;
         if(params.hit_ball) {
@@ -547,20 +556,16 @@ class HPolytope {
         NT* sum_nom_data = sum_nom.data();
         NT* Av_data = Av.data();
 
-        for (int i = 0; i < m; i++) {
-            if (*Av_data == NT(0)) {
-                //std::cout<<"div0"<<std::endl;
-                ;
-            } else {
-                lamda = *sum_nom_data / *Av_data;
-                if (lamda < min_plus && lamda > 0) {
-                    min_plus = lamda;
-                    facet = i;
-                    params.inner_vi_ak = *Av_data;
-                }
+        for (int i = 0; i < m; ++i, ++Av_data, ++sum_nom_data){
+            if (i == skip) continue;    
+            if (*Av_data == NT(0)) continue;
+
+            lamda = *sum_nom_data / *Av_data;
+            if (lamda < min_plus && lamda > 0) {
+                min_plus = lamda;
+                facet = i;
+                params.inner_vi_ak = *Av_data;
             }
-            Av_data++;
-            sum_nom_data++;
         }
         params.facet_prev = facet;
         return std::pair<NT, int>(min_plus, facet);
@@ -641,6 +646,7 @@ class HPolytope {
         NT lamda = 0;
         VT sum_nom;
         int m = num_of_hyperplanes(), facet;
+        int skip = params.facet_prev;
 
         Ar.noalias() += lambda_prev*Av;
         sum_nom.noalias() = b - Ar;
@@ -649,24 +655,20 @@ class HPolytope {
         NT* sum_nom_data = sum_nom.data();
         NT* Av_data = Av.data();
 
-        for (int i = 0; i < m; i++) {
-            if (*Av_data == NT(0)) {
-                //std::cout<<"div0"<<std::endl;
-                ;
-            } else {
-                lamda = *sum_nom_data / *Av_data;
-                if (lamda < min_plus && lamda > 0) {
-                    min_plus = lamda;
-                    facet = i;
-                    params.inner_vi_ak = *Av_data;
-                }
+        for (int i = 0; i < m; ++i, ++Av_data, ++sum_nom_data) {
+            if (i == skip) continue;
+            if (*Av_data == NT(0)) continue;
+            lamda = *sum_nom_data / *Av_data;
+            if (lamda < min_plus && lamda > 0) {
+                min_plus = lamda;
+                facet = i;
+                params.inner_vi_ak = *Av_data;
             }
-            Av_data++;
-            sum_nom_data++;
         }
         params.facet_prev = facet;
         return std::pair<NT, int>(min_plus, facet);
     }
+
 
     //-----------------------------------------------------------------------------------//
 
@@ -1086,4 +1088,4 @@ class HPolytope {
     }
 };
 
-#endif
+#endif

include/diagnostics/scaling_ratio.hpp

@@ -0,0 +1,126 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2025 Vissarion Fisikopoulos
+// Copyright (c) 2018-2025 Apostolos Chalkis
+// Copyright (c) 2025-2025 Iva Janković
+
+// Contributed and/or modified by Iva Janković, as part of Google Summer of Code 2025 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+#ifndef DIAGNOSTICS_SCALING_RATIO_HPP
+#define DIAGNOSTICS_SCALING_RATIO_HPP
+
+#include <tuple>
+#include <vector>
+#include <cmath>
+#include <cstdlib>
+
+template<typename Polytope>
+std::tuple<typename Polytope::VT, typename Polytope::MT,typename Polytope::VT, typename Polytope::VT>
+scaling_ratio_boundary_test(const Polytope&  P,
+                            const typename Polytope::MT& samples,
+                            const typename Polytope::NT tol = 1e-10,
+                            const typename Polytope::NT min_ratio = 0.01) 
+{
+    using VT = typename Polytope::VT;
+    using MT = typename Polytope::MT;
+    using NT = typename Polytope::NT;
+
+    const int dim = P.dimension();
+    const int m = P.num_of_hyperplanes();
+    const unsigned n_samp = static_cast<unsigned>(samples.cols());
+
+    VT scale(10); //we scale 10 times
+    MT coverage(m, 10);
+    std::vector<int> facet_id(n_samp, -1);
+    const auto A_full = P.get_mat();
+    const auto b_full = P.get_vec();
+
+    for (int i = 0; i < n_samp; ++i) {
+        auto Aq = A_full * samples.col(i);
+
+        for (size_t k = 0; k < m; ++k) {
+            if (std::abs(Aq[k] - b_full[k]) < tol) 
+            {
+                facet_id[i] = static_cast<int>(k);
+                break;
+            }
+        }
+    }
+
+    for (int f = 0; f < m; ++f) 
+    {
+        // Samples S on facet f
+        std::vector<int> S;
+        S.reserve(n_samp / m);
+        for (unsigned i = 0; i < n_samp; ++i) {
+            if (facet_id[i] == f)
+                S.push_back(static_cast<int>(i));
+        }
+
+        const double ratio = static_cast<double>(S.size()) / n_samp;
+        if (ratio < min_ratio)  continue;
+
+        // Finding the center
+        VT p = VT::Zero(dim);
+        for (int idx : S) p += samples.col(idx);
+        p /= static_cast<double>(S.size());
+
+        // Looping over scale factor
+        for (int k = 0; k < 10; ++k) 
+        {
+            NT step = 0.1 * (k+1);
+            NT x = std::pow(step, 1.0 / dim);
+            // Local copy of polytope for each scaling
+            Polytope P_loc = P;
+
+            //Shifting and scaling
+            P_loc.shift(p);
+            MT T = (1.0 / x) * MT::Identity(dim, dim);
+            P_loc.linear_transformIt(T);
+
+            //Parameters of new polytope
+            const auto& A_sh = P_loc.get_mat();
+            const auto& b_sh = P_loc.get_vec();
+
+            // Points still in facet
+            unsigned survivors = 0;
+            for (int idx : S) {
+
+                const VT q_shift = samples.col(idx) - p;
+                bool inside = true;
+
+                for (int j = 0; j < A_sh.rows(); ++j) {
+                    if (j == f) continue;
+                    if (A_sh.row(j).dot(q_shift) - b_sh[j] > tol) {
+                        inside = false; break;
+                    }
+                }
+                if (inside) ++survivors;
+            }
+
+            coverage(f, k) = double(survivors) / double(S.size());
+            scale[k]=std::pow(x, dim);
+        }
+    }
+
+    int K = scale.size();           // number of scaling factors
+    VT max_dev(m), avg_dev(m);      // Vectors of average and maximum deviations
+
+    for (int f = 0; f < m; ++f) {
+        double sumd = 0.0;
+        double maxd = 0.0;
+        for (int k = 0; k < K; ++k) {
+            double d = std::abs(coverage(f, k) - scale[k])* 100.0; // in percentage
+            sumd += d;
+            if (d > maxd) maxd = d;
+        }
+        avg_dev[f] = sumd/K;
+        max_dev[f] = maxd;      
+    }
+
+    return {scale, coverage, max_dev, avg_dev};
+}
+
+#endif // DIAGNOSTICS_SCALING_RATIO_HPP

include/preprocess/feasible_point.hpp

@@ -30,5 +30,41 @@ VT compute_feasible_point(MT const& A, VT const& b)
     return x;
 }
 
+template <typename Point, typename Polytope, typename RNG>
+std::pair<typename Polytope::VT,int> compute_feasible_boundary_point(Polytope const& P,
+                                                                     RNG& rng,
+                                                                     typename Point::FT eps)
+{
+    using VT = typename Polytope::VT;
+    using NT = typename Point::FT;
+    const std::size_t m = P.num_of_hyperplanes();
+
+    // Find the interior point
+    VT r = compute_feasible_point(P.get_mat(), P.get_vec());
+
+    // Random ray 
+    const int dim = P.dimension();
+    Point v_pt = GetDirection<Point>::apply(dim, rng);
+    VT v = v_pt.getCoefficients();
+
+    // First‐hit oracle
+    VT Ar(m), Av(m);
+    struct Params { NT inner_vi_ak; int facet_prev; };
+    Params params;
+
+    auto [lambda_min, facet] = P.line_first_positive_intersect(r, v, Ar, Av, params);
+
+    // Checks
+    if (!std::isfinite(lambda_min) || lambda_min <= eps || facet < 0)
+        throw std::runtime_error("Failed to hit boundary!!!");
+
+    // Compute boundry point + final check
+    VT x = r + lambda_min * v;
+    if ((P.get_mat() * x - P.get_vec()).maxCoeff() > eps)
+        throw std::runtime_error("Boundary point violates constraints!!!");
+
+    return std::pair<VT,int>(x, facet);    
+}
+
 
 #endif