Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <iomanip>
- #include <fstream>
- #include <array>
- #include <random>
- #include <cmath>
- #include <cassert>
- #include <limits>
- typedef std::array<double, 2> vector;
- typedef std::array<vector, 6> point_list;
- const auto pi = 4*std::atan(1);
- double random_uniform(double min, double max)
- {
- static std::default_random_engine generator;
- static std::uniform_real_distribution<double> unif_dist; // min = 0, max = 1
- return min + (max - min)*unif_dist(generator);
- }
- double distance(const vector& a, const vector& b)
- {
- return std::sqrt(std::pow(a[0] - b[0], 2) + std::pow(a[1] - b[1], 2));
- }
- bool inside_semicircle(const vector& p)
- {
- return p[1] >= 0 && (distance(p, {0.0, 0.0}) <= 1);
- }
- // raondomly move a single point
- void jostle_point(vector& p, double jostle_distance)
- {
- while(true)
- {
- auto distance = random_uniform(0, jostle_distance);
- auto angle = random_uniform(0, 2*pi);
- auto dx = distance*std::cos(angle);
- auto dy = distance*std::sin(angle);
- auto np = p;
- np[0] += dx;
- np[1] += dy;
- if(inside_semicircle(np))
- {
- p = np;
- break;
- }
- }
- }
- // Move all points randomly and independantly
- void jostle_points_free(point_list& pl, double jostle_distance)
- {
- for(auto& p : pl)
- {
- jostle_point(p, jostle_distance);
- }
- }
- // Write out point list to output stream
- void write_list(const point_list& pl, std::ostream& os)
- {
- os << "p = [\n";
- for(const auto& p : pl)
- {
- for(auto e : p)
- {
- os << std::setprecision(std::numeric_limits<double>::max_digits10) << e << ' ';
- }
- os << ";\n";
- }
- os << "];\n=====================================================" << std::endl;
- }
- // write point list to console and file
- std::ofstream ofs("result.txt");
- void write_list(const point_list& pl)
- {
- write_list(pl, std::cout);
- write_list(pl, ofs);
- }
- // Given a radius associated with points i, j, k
- // in the point list, check if any othe points
- // fall within a radius of the given point (usually
- // the circumcenter of i, j, k).
- bool is_delaunay(const point_list& pl,
- const vector& point,
- double radius,
- size_t i, size_t j, size_t k = -1)
- {
- if( ! std::isfinite(radius))
- {
- return false;
- }
- for(size_t n = 0; n < pl.size(); ++n)
- {
- if((n == i) || (n == j) || (n == k))
- {
- continue;
- }
- if(distance(point, pl[n]) < radius)
- {
- return false;
- }
- }
- return true;
- }
- vector midpoint(const vector& a, const vector& b)
- {
- return {(a[0] + b[0])/2.0, (a[1] + b[1])/2.0};
- }
- // return a vector perpendicular to the line joining a and b
- vector perpendicular(const vector& a, const vector& b)
- {
- return {a[1] - b[1], b[0] - a[0]};
- }
- // dot product
- double dot(const vector& a, const vector& b)
- {
- return (a[0]*b[0]) + (a[1]*b[1]);
- }
- // Note: all maximum distance functions return -1.0
- // to indicate that an invalid or irrelevant distance
- // was the result. Since a maximum distance is being sought,
- // this value will never be used.
- // local maximum distance from the line y = 0 (ye0)
- double max_distance_to_ye0(size_t i, size_t j, const point_list& pl)
- {
- auto mid = midpoint(pl[i], pl[j]);
- auto perp = perpendicular(pl[i], pl[j]);
- // find point where mid + k*perp intersects y = 0
- auto t = -mid[1]/perp[1];
- auto x = mid[0] + t*perp[0];
- auto line_intersect = vector{x, 0.0};
- auto radius = distance(line_intersect, pl[i]);
- if(is_delaunay(pl, line_intersect, radius, i, j) &&
- inside_semicircle(line_intersect))
- {
- return radius;
- }
- return -1.0;
- }
- // Local maximum distance to the circle (r = 1)
- double max_distance_to_re1(size_t i, size_t j, const point_list& pl)
- {
- auto mid = midpoint(pl[i], pl[j]);
- auto perp = perpendicular(pl[i], pl[j]);
- auto dmp = dot(mid, perp);
- double dpp = dot(perp, perp);
- double dmm = dot(mid, mid);
- auto det = std::sqrt(dmp*dmp - dpp*(dmm-1.0));
- auto k1 = (-dmp + det)/dpp;
- auto k2 = (-dmp - det)/dpp;
- vector circle_intersect1 = {mid[0] + k1*perp[0], mid[1] + k1*perp[1]};
- vector circle_intersect2 = {mid[0] + k2*perp[0], mid[1] + k2*perp[1]};
- auto radius1 = distance(circle_intersect1, pl[i]);
- auto radius2 = distance(circle_intersect2, pl[i]);
- if(circle_intersect1[1] < 0.0)
- {
- radius1 = 3.0; // larger than any valid radius
- }
- if(circle_intersect2[1] < 0.0)
- {
- radius2 = 3.0;
- } // these two conditions cannot be true at the same time
- auto circle_intersect = (radius1 < radius2 ? circle_intersect1 : circle_intersect2);
- auto radius = (radius1 < radius2 ? radius1 : radius2);
- if( ! is_delaunay(pl, circle_intersect, radius, i, j))
- {
- return -1.0;
- }
- return radius;
- }
- // Find point equidistant from three points
- vector circumcenter(size_t i, size_t j, size_t k, const point_list& pl)
- {
- auto m0 = midpoint(pl[i], pl[j]);
- auto m1 = midpoint(pl[j], pl[k]);
- auto p0 = perpendicular(pl[i], pl[j]);
- auto p1 = perpendicular(pl[j], pl[k]);
- auto t = ((m1[1] - m0[1])*p0[0] + (m0[0] - m1[0])*p0[1])/(p1[0]*p0[1] - p1[1]*p0[0]);
- return {m1[0] + t*p1[0], m1[1] + t*p1[1]};
- }
- // return the maximum distance from three points while still occupying
- // the triangle formed by the three points
- double circumcircle_radius(size_t i, size_t j, size_t k, const point_list& pl)
- {
- auto circumc = circumcenter(i, j, k, pl);
- auto radius = distance(pl[i], circumc);
- assert(distance(pl[i], circumc) - distance(pl[j], circumc) < 1e-4);
- assert(distance(pl[k], circumc) - distance(pl[j], circumc) < 1e-4);
- if(std::isfinite(radius) &&
- is_delaunay(pl, circumc, radius, i, j, k) &&
- inside_semicircle(circumc))
- {
- return radius;
- }
- return -1.0;
- }
- // find the maximum distance to any point in the point list
- // using the above maximum distance functions.
- double maximum_minimum_distance(const point_list& pl)
- {
- auto maximum_distance = 0.0;
- // Find maximum distance inside convex hull of points
- for(size_t i = 0; i < pl.size(); ++i)
- {
- for(size_t j = i + 1; j < pl.size(); ++j)
- {
- for(size_t k = j + 1; k < pl.size(); ++k)
- {
- auto del_radius = circumcircle_radius(i, j, k, pl);
- if(del_radius > maximum_distance)
- {
- maximum_distance = del_radius;
- }
- }
- }
- }
- // For each pair of points, find the maximum distance to the
- // semicircular border and straight edge
- for(size_t i = 0; i < pl.size(); ++i)
- {
- for(size_t j = i + 1; j < pl.size(); ++j)
- {
- // find point where mid + k*perp intersects y = 0
- auto max_distance_to_line = max_distance_to_ye0(i, j, pl);
- if(max_distance_to_line > maximum_distance)
- {
- maximum_distance = max_distance_to_line;
- }
- // find point where mid + k*perp intersects x^2 + y^2 = 1
- auto max_distance_to_circle = max_distance_to_re1(i, j, pl);
- if(max_distance_to_circle > maximum_distance)
- {
- maximum_distance = max_distance_to_circle;
- }
- }
- }
- // check each point for distance to corners (-1, 0) and (1, 0)
- auto minimum_distance_to_left_corner = 2.0;
- for(const auto& p : pl)
- {
- auto d = distance(p, {-1.0, 0.0});
- if(d < minimum_distance_to_left_corner)
- {
- minimum_distance_to_left_corner = d;
- }
- }
- if(minimum_distance_to_left_corner > maximum_distance)
- {
- maximum_distance = minimum_distance_to_left_corner;
- }
- auto minimum_distance_to_right_corner = 2.0;
- for(const auto& p : pl)
- {
- auto d = distance(p, {1.0, 0.0});
- if(d < minimum_distance_to_right_corner)
- {
- minimum_distance_to_right_corner = d;
- }
- }
- if(minimum_distance_to_right_corner > maximum_distance)
- {
- maximum_distance = minimum_distance_to_right_corner;
- }
- return maximum_distance;
- }
- int main()
- {
- point_list pl;
- for(size_t i = 0; i < pl.size(); ++i)
- {
- pl[i] = {0.0, 0.5};
- }
- auto jostle_distance = 2.0;
- auto min_max_min_dis = 3.0;
- auto fail_count = 0;
- auto fail_threshold = 1000000;
- auto success = false;
- while(jostle_distance < 4.0)
- {
- auto new_pl = pl;
- jostle_points_free(new_pl, jostle_distance);
- double max_min_dis = maximum_minimum_distance(new_pl);
- if(max_min_dis < min_max_min_dis)
- {
- pl = new_pl;
- min_max_min_dis = max_min_dis;
- fail_count = 0;
- success = true;
- std::cout << "r = " << std::setprecision(std::numeric_limits<double>::max_digits10) << max_min_dis
- << " % Jostle = " << jostle_distance << std::endl;
- ofs << "r = " << std::setprecision(std::numeric_limits<double>::max_digits10) << max_min_dis
- << " % Jostle = " << jostle_distance << std::endl;
- write_list(pl);
- }
- else
- {
- if(++fail_count >= fail_threshold)
- {
- if(success)
- {
- // If new low found, take smaller steps
- jostle_distance *= 0.5;
- }
- else
- {
- // if no new low found, take larger steps to escape local minimum
- jostle_distance *= 2.0;
- }
- // don't let the jostle_distance get so small that adding
- // it to a position doesn't change anything
- while(jostle_distance < std::numeric_limits<double>::epsilon())
- {
- jostle_distance *= 20.0;
- }
- fail_count = 0;
- success = false;
- std::cout << "New jostle distance = " << jostle_distance << std::endl;
- }
- }
- }
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement