SHARE
TWEET

semicircle points

a guest Sep 23rd, 2015 132 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include <iomanip>
  3. #include <fstream>
  4. #include <array>
  5. #include <random>
  6. #include <cmath>
  7. #include <cassert>
  8. #include <limits>
  9.  
  10. typedef std::array<double, 2> vector;
  11. typedef std::array<vector, 6> point_list;
  12.  
  13. const auto pi = 4*std::atan(1);
  14.  
  15. double random_uniform(double min, double max)
  16. {
  17.     static std::default_random_engine generator;
  18.     static std::uniform_real_distribution<double> unif_dist; // min = 0, max = 1
  19.  
  20.     return min + (max - min)*unif_dist(generator);
  21. }
  22.  
  23. double distance(const vector& a, const vector& b)
  24. {
  25.     return std::sqrt(std::pow(a[0] - b[0], 2) + std::pow(a[1] - b[1], 2));
  26. }
  27.  
  28. bool inside_semicircle(const vector& p)
  29. {
  30.     return p[1] >= 0 && (distance(p, {0.0, 0.0}) <= 1);
  31. }
  32.  
  33. // randomly moves points in point_list while obeying symmetry constraints
  34. void jostle_points_constrained(point_list& pl, double jostle_distance)
  35. {
  36.     size_t index;
  37.     if(pl.size() == 6)
  38.     {
  39.         index = 2;
  40.     }
  41.     else
  42.     {
  43.         assert(false);
  44.     }
  45.  
  46.     // First point(s) stay on the y-axis
  47.     for(auto i = 0; i < index; ++i)
  48.     {
  49.         while(true)
  50.         {
  51.             auto distance = random_uniform(-jostle_distance, jostle_distance);
  52.             auto np = pl[i];
  53.             np[0] = 0.0;
  54.             np[1] += distance;
  55.             if(inside_semicircle(np))
  56.             {
  57.                 pl[i] = np;
  58.                 break;
  59.             }
  60.         }
  61.     }
  62.  
  63.     // Other points are mirrored across the y-axis
  64.     for(size_t i = index; i < pl.size(); i += 2)
  65.     {
  66.         while(true)
  67.         {
  68.             auto distance = random_uniform(0, jostle_distance);
  69.             auto angle = random_uniform(0, 2*pi);
  70.  
  71.             auto dx = distance*std::cos(angle);
  72.             auto dy = distance*std::sin(angle);
  73.  
  74.             auto np = pl[i];
  75.             np[0] += dx;
  76.             np[1] += dy;
  77.             if(inside_semicircle(np))
  78.             {
  79.                 pl[i] = np;
  80.                 pl[i+1] = np;
  81.                 pl[i+1][0] *= -1.0; // reflect in y-axis
  82.                 break;
  83.             }
  84.         }
  85.     }
  86. }
  87.  
  88. // raondomly move a single point
  89. void jostle_point(vector& p, double jostle_distance)
  90. {
  91.     while(true)
  92.     {
  93.         auto distance = random_uniform(0, jostle_distance);
  94.         auto angle = random_uniform(0, 2*pi);
  95.  
  96.         auto dx = distance*std::cos(angle);
  97.         auto dy = distance*std::sin(angle);
  98.  
  99.         auto np = p;
  100.         np[0] += dx;
  101.         np[1] += dy;
  102.         if(inside_semicircle(np))
  103.         {
  104.             p = np;
  105.             break;
  106.         }
  107.     }
  108. }
  109.  
  110. // Move all points randomly and independantly
  111. void jostle_points_free(point_list& pl, double jostle_distance)
  112. {
  113.     for(auto& p : pl)
  114.     {
  115.         jostle_point(p, jostle_distance);
  116.     }
  117. }
  118.  
  119. // Write out point list to output stream
  120. void write_list(const point_list& pl, std::ostream& os)
  121. {
  122.     os << "p = [\n";
  123.     for(const auto& p : pl)
  124.     {
  125.         for(auto e : p)
  126.         {
  127.             os << std::setprecision(std::numeric_limits<double>::max_digits10) << e << ' ';
  128.         }
  129.         os << ";\n";
  130.     }
  131.     os << "];\n=====================================================" << std::endl;
  132. }
  133.  
  134. // write point list to console and file
  135. std::ofstream ofs("result.txt");
  136. void write_list(const point_list& pl)
  137. {
  138.     write_list(pl, std::cout);
  139.     write_list(pl, ofs);
  140. }
  141.  
  142. // Given a radius associated with points i, j, k
  143. // in the point list, check if any othe points
  144. // fall within a radius of the given point (usually
  145. // the circumcenter of i, j, k).
  146. bool is_delaunay(const point_list& pl,
  147.                  const vector& point,
  148.                  double radius,
  149.                  size_t i, size_t j, size_t k = -1)
  150. {
  151.     if( ! std::isfinite(radius))
  152.     {
  153.         return false;
  154.     }
  155.  
  156.     for(size_t n = 0; n < pl.size(); ++n)
  157.     {
  158.         if((n == i) || (n == j) || (n == k))
  159.         {
  160.             continue;
  161.         }
  162.  
  163.         if(distance(point, pl[n]) < radius)
  164.         {
  165.             return false;
  166.         }
  167.     }
  168.  
  169.     return true;
  170. }
  171.  
  172. vector midpoint(const vector& a, const vector& b)
  173. {
  174.     return {(a[0] + b[0])/2.0, (a[1] + b[1])/2.0};
  175. }
  176.  
  177. // return a vector perpendicular to the line joining a and b
  178. vector perpendicular(const vector& a, const vector& b)
  179. {
  180.     return {a[1] - b[1], b[0] - a[0]};
  181. }
  182.  
  183. // dot product
  184. double dot(const vector& a, const vector& b)
  185. {
  186.     return (a[0]*b[0]) + (a[1]*b[1]);
  187. }
  188.  
  189. // Note: all maximum distance functions return -1.0
  190. // to indicate that an invalid or irrelevant distance
  191. // was the result. Since a maximum distance is being sought,
  192. // this value will never be used.
  193.  
  194. // local maximum distance from the line y = 0 (ye0)
  195. double max_distance_to_ye0(size_t i, size_t j, const point_list& pl)
  196. {
  197.     auto mid = midpoint(pl[i], pl[j]);
  198.     auto perp = perpendicular(pl[i], pl[j]);
  199.  
  200.     // find point where mid + k*perp intersects y = 0
  201.     auto t = -mid[1]/perp[1];
  202.     auto x = mid[0] + t*perp[0];
  203.     auto line_intersect = vector{x, 0.0};
  204.     auto radius = distance(line_intersect, pl[i]);
  205.     if(is_delaunay(pl, line_intersect, radius, i, j) &&
  206.        inside_semicircle(line_intersect))
  207.     {
  208.         return radius;
  209.     }
  210.  
  211.     return -1.0;
  212. }
  213.  
  214. // Local maximum distance to the circle (r = 1)
  215. double max_distance_to_re1(size_t i, size_t j, const point_list& pl)
  216. {
  217.     auto mid = midpoint(pl[i], pl[j]);
  218.     auto perp = perpendicular(pl[i], pl[j]);
  219.     auto dmp = dot(mid, perp);
  220.     double dpp = dot(perp, perp);
  221.     double dmm = dot(mid, mid);
  222.     auto det = std::sqrt(dmp*dmp - dpp*(dmm-1.0));
  223.     auto k1 = (-dmp + det)/dpp;
  224.     auto k2 = (-dmp - det)/dpp;
  225.     vector circle_intersect1 = {mid[0] + k1*perp[0], mid[1] + k1*perp[1]};
  226.     vector circle_intersect2 = {mid[0] + k2*perp[0], mid[1] + k2*perp[1]};
  227.     auto radius1 = distance(circle_intersect1, pl[i]);
  228.     auto radius2 = distance(circle_intersect2, pl[i]);
  229.     if(circle_intersect1[1] < 0.0)
  230.     {
  231.         radius1 = 3.0; // larger than any valid radius
  232.     }
  233.     if(circle_intersect2[1] < 0.0)
  234.     {
  235.         radius2 = 3.0;
  236.     } // these two conditions cannot be true at the same time
  237.  
  238.     auto circle_intersect = (radius1 < radius2 ? circle_intersect1 : circle_intersect2);
  239.     auto radius = (radius1 < radius2 ? radius1 : radius2);
  240.  
  241.     if( ! is_delaunay(pl, circle_intersect, radius, i, j))
  242.     {
  243.         return -1.0;
  244.     }
  245.  
  246.     return radius;
  247. }
  248.  
  249. // Find point equidistant from three points
  250. vector circumcenter(size_t i, size_t j, size_t k, const point_list& pl)
  251. {
  252.     auto m0 = midpoint(pl[i], pl[j]);
  253.     auto m1 = midpoint(pl[j], pl[k]);
  254.  
  255.     auto p0 = perpendicular(pl[i], pl[j]);
  256.     auto p1 = perpendicular(pl[j], pl[k]);
  257.  
  258.     auto t = ((m1[1] - m0[1])*p0[0] + (m0[0] - m1[0])*p0[1])/(p1[0]*p0[1] - p1[1]*p0[0]);
  259.  
  260.     return {m1[0] + t*p1[0], m1[1] + t*p1[1]};
  261. }
  262.  
  263. // return the maximum distance from three points while still occupying
  264. // the triangle formed by the three points
  265. double circumcircle_radius(size_t i, size_t j, size_t k, const point_list& pl)
  266. {
  267.     auto circumc = circumcenter(i, j, k, pl);
  268.     auto radius = distance(pl[i], circumc);
  269.  
  270.     assert(distance(pl[i], circumc) - distance(pl[j], circumc) < 1e-4);
  271.     assert(distance(pl[k], circumc) - distance(pl[j], circumc) < 1e-4);
  272.  
  273.     if(std::isfinite(radius) &&
  274.        is_delaunay(pl, circumc, radius, i, j, k) &&
  275.        inside_semicircle(circumc))
  276.     {
  277.         return radius;
  278.     }
  279.  
  280.     return -1.0;
  281. }
  282.  
  283. // find the maximum distance to any point in the point list
  284. // using the above maximum distance functions.
  285. double maximum_minimum_distance(const point_list& pl)
  286. {
  287.     auto maximum_distance = 0.0;
  288.  
  289.     // Find maximum distance inside convex hull of points
  290.     for(size_t i = 0; i < pl.size(); ++i)
  291.     {
  292.         for(size_t j = i + 1; j < pl.size(); ++j)
  293.         {
  294.             for(size_t k = j + 1; k < pl.size(); ++k)
  295.             {
  296.                 auto del_radius = circumcircle_radius(i, j, k, pl);
  297.  
  298.                 if(del_radius > maximum_distance)
  299.                 {
  300.                     maximum_distance = del_radius;
  301.                 }
  302.             }
  303.         }
  304.     }
  305.  
  306.     // For each pair of points, find the maximum distance to the
  307.     // semicircular border and straight edge
  308.     for(size_t i = 0; i < pl.size(); ++i)
  309.     {
  310.         for(size_t j = i + 1; j < pl.size(); ++j)
  311.         {
  312.             // find point where mid + k*perp intersects y = 0
  313.             auto max_distance_to_line = max_distance_to_ye0(i, j, pl);
  314.             if(max_distance_to_line > maximum_distance)
  315.             {
  316.                 maximum_distance = max_distance_to_line;
  317.             }
  318.  
  319.             // find point where mid + k*perp intersects x^2 + y^2 = 1
  320.             auto max_distance_to_circle = max_distance_to_re1(i, j, pl);
  321.             if(max_distance_to_circle > maximum_distance)
  322.             {
  323.                 maximum_distance = max_distance_to_circle;
  324.             }
  325.         }
  326.     }
  327.  
  328.     // check each point for distance to corners (-1, 0) and (1, 0)
  329.     auto minimum_distance_to_left_corner = 2.0;
  330.     for(const auto& p : pl)
  331.     {
  332.         auto d = distance(p, {-1.0, 0.0});
  333.         if(d < minimum_distance_to_left_corner)
  334.         {
  335.             minimum_distance_to_left_corner = d;
  336.         }
  337.     }
  338.     if(minimum_distance_to_left_corner > maximum_distance)
  339.     {
  340.         maximum_distance = minimum_distance_to_left_corner;
  341.     }
  342.  
  343.     auto minimum_distance_to_right_corner = 2.0;
  344.     for(const auto& p : pl)
  345.     {
  346.         auto d = distance(p, {1.0, 0.0});
  347.         if(d < minimum_distance_to_right_corner)
  348.         {
  349.             minimum_distance_to_right_corner = d;
  350.         }
  351.     }
  352.     if(minimum_distance_to_right_corner > maximum_distance)
  353.     {
  354.         maximum_distance = minimum_distance_to_right_corner;
  355.     }
  356.  
  357.     return maximum_distance;
  358. }
  359.  
  360. int main()
  361. {
  362.     point_list pl;
  363.     for(size_t i = 0; i < pl.size(); ++i)
  364.     {
  365.         pl[i] = {0.0, 0.5};
  366.     }
  367.     auto jostle_distance = 1e-4;
  368.     auto min_max_min_dis = 3.0;
  369.     auto fail_count = 0;
  370.     auto fail_threshold = 1000000;
  371.     auto success = false;
  372.  
  373.     while(jostle_distance < 4.0)
  374.     {
  375.         auto new_pl = pl;
  376.         jostle_points_constrained(new_pl, jostle_distance);
  377.         double max_min_dis = maximum_minimum_distance(new_pl);
  378.         if(max_min_dis < min_max_min_dis)
  379.         {
  380.             pl = new_pl;
  381.             min_max_min_dis = max_min_dis;
  382.             fail_count = 0;
  383.             success = true;
  384.  
  385.             std::cout << "r = " << std::setprecision(std::numeric_limits<double>::max_digits10) << max_min_dis
  386.                       << " % Jostle = " << jostle_distance << std::endl;
  387.             ofs       << "r = " << std::setprecision(std::numeric_limits<double>::max_digits10) << max_min_dis
  388.                       << " % Jostle = " << jostle_distance << std::endl;
  389.             write_list(pl);
  390.         }
  391.         else
  392.         {
  393.             if(++fail_count >= fail_threshold)
  394.             {
  395.                 if(success)
  396.                 {
  397.                     // If new low found, take smaller steps
  398.                     jostle_distance *= 0.5;
  399.                 }
  400.                 else
  401.                 {
  402.                     // if no new low found, take larger steps to escape local minimum
  403.                     jostle_distance *= 2.0;
  404.                 }
  405.  
  406.                 // don't let the jostle_distance get so small that adding
  407.                 // it to a position doesn't change anything
  408.                 while(jostle_distance < std::numeric_limits<double>::epsilon())
  409.                 {
  410.                     jostle_distance *= 20.0;
  411.                 }
  412.  
  413.                 fail_count = 0;
  414.                 success = false;
  415.                 std::cout << "New jostle distance = " << jostle_distance << std::endl;
  416.             }
  417.         }
  418.     }
  419.  
  420.     return 0;
  421. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top