semicircle points

#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);
}

// randomly moves points in point_list while obeying symmetry constraints
void jostle_points_constrained(point_list& pl, double jostle_distance)
{
    size_t index;
    if(pl.size() == 6)
    {
        index = 2;
    }
    else
    {
        assert(false);
    }

    // First point(s) stay on the y-axis
    for(auto i = 0; i < index; ++i)
    {
        while(true)
        {
            auto distance = random_uniform(-jostle_distance, jostle_distance);
            auto np = pl[i];
            np[0] = 0.0;
            np[1] += distance;
            if(inside_semicircle(np))
            {
                pl[i] = np;
                break;
            }
        }
    }

    // Other points are mirrored across the y-axis
    for(size_t i = index; i < pl.size(); i += 2)
    {
        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 = pl[i];
            np[0] += dx;
            np[1] += dy;
            if(inside_semicircle(np))
            {
                pl[i] = np;
                pl[i+1] = np;
                pl[i+1][0] *= -1.0; // reflect in y-axis
                break;
            }
        }
    }
}

// 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 = 1e-4;
    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_constrained(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;
}