Aumente kamada_kawai_spring_layout usando topologia derivada de GeographicLib ... falhando

Então, depois de revisar, criei uma versão completa.

Havia o bug que inseria 0como vizinho qualquer nó (o inicializador {{}, {}}cria um conjunto com o único elemento 0, porque ambos {}os valores são inicializados para esse valor).

Além disso, a rede armazena vizinhos em ambas as direções, mas como o gráfico é undirectedde qualquer forma, isso resulta em todas as arestas sendo duplicadas. Se essa não for a intenção, simplesmente substitua vecSpor setSpara as arestas e as duplicatas serão filtradas automaticamente.

Então, a propriedade edge index é necessária, mas não foi preenchida corretamente. Elas foram todas definidas como 1(1.0, então parece que isso era /intencionado/ como o peso).

Nem o peso. Todos foram definidos como 0.

Acho que é ainda melhor não começar com last_delta = 0.

Por fim, provavelmente está tudo bem, mas eu queria ter certeza de que o last_deltaseria realmente mantido, passando a condição "done" por referência.

Por fim, adicionei uma write_graphvizetapa para renderizar a entrada conforme analisada para inspeção:

Ao vivo em Coliru

#include "rapidjson/document.h"
#include "rapidjson/ostreamwrapper.h"
#include "rapidjson/writer.h"
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/Intersect.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graphviz.hpp>
#include <boost/graph/kamada_kawai_spring_layout.hpp>
#include <boost/graph/random_layout.hpp>
#include <fstream>
#include <iostream>
#include <random>

namespace GL = GeographicLib;

template <typename RandomNumberGenerator = std::minstd_rand>
class geo_topology : public boost::convex_topology<2> {
  public:
    using point_type            = typename convex_topology<2>::point_type;
    using point_difference_type = typename convex_topology<2>::point_difference_type;
    using rand_t                = boost::uniform_01<RandomNumberGenerator, double>;

  private:
    RandomNumberGenerator gen_;
    mutable rand_t        rand_{gen_};
    GL::Geodesic          geodesic_;

  public:
    geo_topology(geo_topology const&)            = delete;
    geo_topology& operator=(geo_topology const&) = delete;

    explicit geo_topology(GL::Geodesic const& g) : geodesic_(g) {}
    geo_topology(RandomNumberGenerator gen) : gen_(std::move(gen)), geodesic_(GL::Geodesic::WGS84()) {}

    point_type random_point() const {
        point_type p;
        p[0] = fmod(rand_(), 180.0) - 90.0;  // latitude.
        p[1] = fmod(rand_(), 360.0) - 180.0; // longitude.
        if (std::isnan(p[0]) or std::isnan(p[1])) {
            std::cerr << "nan in random_point" << std::endl;
        }
        return p;
    }

    point_type bound(point_type a) const {
        a[0] = fmod(a[0] + 90.0, 180.0) - 90.0;   // latitude.
        a[1] = fmod(a[1] + 180.0, 360.0) - 180.0; // longitude.
        if (std::isnan(a[0]) or std::isnan(a[1])) {
            std::cerr << "nan in bound" << std::endl;
        }
        return a;
    }

    double distance_from_boundary(point_type /*a*/) const {
        return 1000.0; // really don't know if this means anything.
    }

    point_type center() const {
        point_type result;
        result[0] = 90.0; // latitude.
        result[1] = 0.0;
        return result;
    }

    point_type origin() const {
        point_type result;
        result[0] = 0.0; // latitude.
        result[1] = 0.0;
        return result;
    }

    point_difference_type extent() const {
        point_difference_type result;
        result[0] = geodesic_.EquatorialRadius(); // dist.
        result[1] = 360.0;                        // r_fwd_azi
        return result;
    }

    double distance(point a, point b) const {
        double dist = 0.0;
        double r_fwd_azi, r_rev_azi;
        geodesic_.Inverse(a[0], a[1], b[0], b[1], dist, r_fwd_azi, r_rev_azi);
        if (std::isnan(dist)) {
            std::cerr << "dist shows nan distance" << std::endl;
        }
        return dist;
    }

    point move_position_toward(point a, double fraction, point b) const {
        point  c;
        double dist = 0.0;
        double r_fwd_azi, r_rev_azi;
        geodesic_.Inverse(a[0], a[1], b[0], b[1], dist, r_fwd_azi, r_rev_azi);
        GL::GeodesicLine gl = geodesic_.Line(a[0], a[1], r_fwd_azi);
        gl.Position(fraction * dist, c[0], c[1]);
        if (std::isnan(c[0]) or std::isnan(c[1])) {
            std::cerr << "nan in move_position_toward" << std::endl;
        }
        return c;
    }

    point_difference difference(point a, point b) const { // a-b
        point_difference result;
        double           r_rev_azi;
        geodesic_.Inverse(a[0], a[1], b[0], b[1], result[0], result[1], r_rev_azi);
        if (std::isnan(result[0]) or std::isnan(result[1])) {
            std::cerr << "nan in point_difference" << std::endl;
        }
        return result;
    }

    point adjust(point a, point_difference delta) const {
        point result; // point_difference = dist/fwd_azi
        if (std::isnan(delta[0]) or std::isnan(delta[1])) {
            //          std::cerr << "nan in adjust.delta for point at " << a[0] << "," << a[1] << std::endl;
            point b = random_point();
            result  = move_position_toward(a, 0.000001, b);
        } else {
            GL::GeodesicLine gl = geodesic_.Line(a[0], a[1], delta[1]);
            gl.Position(delta[0], result[0], result[1]);
            if (std::isnan(result[0]) or std::isnan(result[1])) {
                std::cerr << "nan in adjust" << std::endl;
            }
        }
        return result;
    }

    point pointwise_min(point a, point b) const {
        point  result;
        double a_dist, b_dist = 0.0;
        double r_fwd_azi, r_rev_azi;
        geodesic_.Inverse(90.0, 0.0, a[0], a[1], a_dist, r_fwd_azi, r_rev_azi);
        geodesic_.Inverse(90.0, 0.0, b[0], b[1], b_dist, r_fwd_azi, r_rev_azi);
        result = a_dist <= b_dist ? a : b;
        if (std::isnan(result[0]) or std::isnan(result[1])) {
            std::cerr << "nan in pointwise_min" << std::endl;
        }
        return result;
    }

    point pointwise_max(point a, point b) const {
        point  result;
        double a_dist, b_dist = 0.0;
        double r_fwd_azi, r_rev_azi;
        geodesic_.Inverse(90.0, 0.0, a[0], a[1], a_dist, r_fwd_azi, r_rev_azi);
        geodesic_.Inverse(90.0, 0.0, b[0], b[1], b_dist, r_fwd_azi, r_rev_azi);
        result = a_dist >= b_dist ? a : b;
        if (std::isnan(result[0]) or std::isnan(result[1])) {
            std::cerr << "nan in pointwise_max" << std::endl;
        }
        return result;
    }

    double norm(point_difference delta) const { return delta[0] / geodesic_.EquatorialRadius(); }

    double volume(point_difference delta) const { return delta[1]; }
};

namespace boost {
    // for graphviz output of node position
    static std::ostream& operator<<(std::ostream& os, geo_topology<>::point const& p) { return os << "(" << p[0] << ", " << p[1] << ")"; }
    static std::istream& operator>>(std::istream& is, geo_topology<>::point&) { is.setstate(std::ios::failbit); return is; } // NOT IMPLEMENTED
} // namespace boost

struct ExampleSpherical {

    static constexpr double r_mux  = std::numbers::pi / 180.0;
    static constexpr double d_mux  = 180.0 / std::numbers::pi;
    static constexpr double radius = 6371008.7714;
    static constexpr double q_circ = radius * std::numbers::pi * 0.5;

    using Topology           = geo_topology<>;
    using Point              = Topology::point_type;
    using EdgeWeightProperty = boost::property<boost::edge_weight_t, double>;

    using Graph = boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS, boost::no_property,
                                        boost::property<boost::edge_index_t, int, EdgeWeightProperty>>;
    using V = boost::graph_traits<Graph>::vertex_descriptor;
    using E = boost::graph_traits<Graph>::edge_descriptor;

    size_t const depth;
    double       line_length;
    using PosVector   = std::vector<Point>;
    using vertex_id_t = unsigned long; // std::size_t;

    std::map<vertex_id_t, std::set<vertex_id_t>> network_;
    std::map<vertex_id_t, Point>                 initial_positions_;
    GL::Geodesic                                 geo;
    Topology                                     topo;

    ExampleSpherical(size_t d) : depth(d), geo(GL::Geodesic::WGS84()), topo(geo) {
        line_length = q_circ / pow(3.0, double(depth + 1));
        std::cout << "depth:" << depth << " line_length:" << line_length << std::endl;
    }

    void load_positions(std::string const filename) {
        //  Load latitude/longitude positions.
        auto                start = std::chrono::system_clock::now();
        rapidjson::Document jpos;
        std::ifstream       pos_file(filename); // dict of node: pos
        std::string         pos_str(std::istreambuf_iterator<char>(pos_file), {});
        jpos.Parse(pos_str.c_str());
        if (!jpos.HasParseError()) {
            for (auto p_it = jpos.MemberBegin(), n_end = jpos.MemberEnd(); p_it != n_end; ++p_it) {
                auto  node = std::stoul(p_it->name.GetString());
                Point p;
                p[0] = p_it->value[0].GetDouble();
                p[1] = p_it->value[1].GetDouble();
                initial_positions_.emplace(node, std::move(p));
            }
        }
        for (auto& [a_id, a_pos] : initial_positions_) {
            for (auto& [b_id, b_pos] : initial_positions_) {
                if (a_id != b_id) {
                    auto dx = topo.difference(a_pos, b_pos);
                    if (dx[0] < line_length) {
                        line_length = dx[0];
                    }
                    if (dx[0] < 0.000001) {
                        std::cout << "found overlapping nodes" << std::endl;
                        exit(0);
                    }
                }
            }
        }
        std::cout << "Minimum Node Distance:" << line_length << std::endl;
        auto end = std::chrono::system_clock::now();
        auto ms  = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
        std::cout << "position load took: " << ms << "ms" << std::endl;
    }

    void load_network(std::string const filename) {
        auto                start = std::chrono::system_clock::now();
        rapidjson::Document jnet;
        std::ifstream net_file(filename); // dict with list "nodes" = each= eg {"id"=0}, and list links=eg
                                          // {"source": 0, "target":2}
        std::string net_str(std::istreambuf_iterator<char>(net_file), {});
        jnet.Parse(net_str.c_str());
        if (!jnet.HasParseError()) {
            for (auto&& n : jnet["nodes"].GetArray())
                network_.insert({n.GetObject()["id"].GetUint64(), {}});

            for (auto&& l : jnet["links"].GetArray()) {
                auto ia = l.GetObject()["source"].GetUint64();
                auto ib = l.GetObject()["target"].GetUint64();
                network_[ia].insert(ib);
                network_[ib].insert(ia);
            }
        }
        auto end = std::chrono::system_clock::now();
        auto ms  = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
        std::cout << "network load took: " << ms << "ms\n";
    }

    void validate_graph() const {
        assert(network_.size() == initial_positions_.size());

        std::set<vertex_id_t> error_nodes;
        for (auto& [id, dd] : network_) {
            if (!initial_positions_.contains(id))
                error_nodes.insert(id); // no position
            if (dd.empty())    // no edges
                error_nodes.insert(id);
            for (auto& q : dd)
                if (!network_.contains(q)) // missing neighbour
                    error_nodes.insert(id);
        }

        if (!error_nodes.empty())
            throw std::runtime_error("files did not reconcile correctly");
    }

    void load_and_validate() {
        load_network("graph_" + std::to_string(depth) + ".json");
        load_positions("graph_ll_pos_" + std::to_string(depth) + ".json");
        validate_graph();
    }

    void save_positions(PosVector& final_pos) const {
        auto                start = std::chrono::system_clock::now();
        rapidjson::Document jpos;
        jpos.SetObject();
        auto& alloc = jpos.GetAllocator();
        for (auto& [id, _] : initial_positions_) {
            rapidjson::Value array(rapidjson::kArrayType);
            array.PushBack(final_pos[id][0], alloc);
            array.PushBack(final_pos[id][1], alloc);
            auto j_id = rapidjson::Value(std::to_string(id).c_str(), alloc);
            jpos.AddMember(j_id, array, alloc);
        }
        std::ofstream ofs{"result_ll_pos_" + std::to_string(depth) + ".json"};
        if (!ofs.is_open()) {
            std::cerr << "Could not open file for writing!\n";
        }
        rapidjson::OStreamWrapper                    osw{ofs};
        rapidjson::Writer<rapidjson::OStreamWrapper> writer2{osw};
        jpos.Accept(writer2);

        auto end = std::chrono::system_clock::now();
        auto ms  = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
        std::cout << "position save took: " << ms << "ms\n";
    }

    struct kamada_kawai_done {
        bool operator()(double delta_p, V /*p*/, Graph const& /*g*/, bool global) {
            std::cout << "kamada_kawai_done " << (global ? "(global) " : "") << delta_p << std::endl;
            if (global) {
                using std::abs;
                auto diff = abs(last_delta - delta_p);
                std::cout << "delta_p: " << delta_p << " diff: " << diff << std::endl;
                last_delta = delta_p;
                return diff < 0.01;
            } else {
                return delta_p < 0.01;
            }
        }
        double last_delta = 0;
    };

    void evaluate() const {
        // std::vector<unsigned long> vertices_(network_.size());
        // for (auto [id, _] : network_)
        //     vertices_.push_back(id);

        Graph     g(network_.size());
        PosVector position(num_vertices(g));
        auto      wmap = get(boost::edge_weight, g);
        auto      eidx = get(boost::edge_index, g);

        for (size_t edge_idx = 0; auto [s, dd] : network_) {
            for (auto d : dd) {
                assert(d != s);
                assert(d < num_vertices(g));
                assert(s < num_vertices(g));
                auto [e, not_duplicate] = add_edge(s, d, g);

                if (not_duplicate) {
                    eidx[e] = edge_idx++; // required to map edges to [0, num_edges(g))
                    wmap[e] = 1.0;        // weight bust all be non-negative
                }
            }
            position[s] = initial_positions_.at(s); // throws if not found.
        }
        auto position_map = make_iterator_property_map(position.begin(), get(boost::vertex_index, g));

        {
            // save graph-viz
            boost::dynamic_properties dp;
            dp.property("node_id", get(boost::vertex_index, g));
            dp.property("pos", position_map);
            dp.property("edge_id", eidx);
            dp.property("weight", wmap);

            std::ofstream dot("graph_" + std::to_string(depth) + ".dot");
            write_graphviz_dp(dot, g, dp);
            // create PNG:
            ::system(("dot -Tpng 'graph_" + std::to_string(depth) + ".dot' -o 'graph_" + std::to_string(depth) + ".png'").c_str());
        }

        kamada_kawai_done condition{line_length / 2.0};
        kamada_kawai_spring_layout(g, position_map, wmap, topo, boost::edge_length(10.0),
                                   std::ref(condition));
        save_positions(position);
    }
};

int main() {
    ExampleSpherical p(0);
    p.load_and_validate();
    p.evaluate();
}

Isto imprime, por exemplo:

depth:0 line_length:3.33585e+06
network load took: 2ms
Minimum Node Distance:3.32011e+06
position load took: 1ms
kamada_kawai_done (global) 0
delta_p: 0 diff: 1.66006e+06
kamada_kawai_done -141.912
kamada_kawai_done (global) -141.9
delta_p: -141.9 diff: 141.9
kamada_kawai_done -141.9
kamada_kawai_done (global) -141.9
delta_p: -141.9 diff: 0
position save took: 0ms

O gráfico de pontos:

graph G {
0  [pos="(0,        180)"];  1  [pos="(0,        -150)"]; 2  [pos="(-30,     180)"];  3    [pos="(30,       180)"];    4  [pos="(0,       150)"];
5  [pos="(-35.2644, -135)"]; 6  [pos="(35.2644,  -135)"]; 7  [pos="(0,       -120)"]; 8    [pos="(-35.2644, 135)"];    9  [pos="(35.2644, 135)"];
10 [pos="(-60,      180)"];  11 [pos="(60,       180)"];  12 [pos="(0,       120)"];  13   [pos="(0,        -90)"];    14 [pos="(-30,     -90)"];
15 [pos="(30,       -90)"];  16 [pos="(-60,      -90)"];  17 [pos="(60,      -90)"];  18   [pos="(-90,      0)"];      19 [pos="(90,      0)"];
20 [pos="(-60,      90)"];   21 [pos="(60,       90)"];   22 [pos="(-30,     90)"];   23   [pos="(30,       90)"];     24 [pos="(0,       90)"];
25 [pos="(0,        -60)"];  26 [pos="(-60,      0)"];    27 [pos="(60,      0)"];    28   [pos="(-35.2644, -45)"];    29 [pos="(35.2644, -45)"];
30 [pos="(0,        60)"];   31 [pos="(-35.2644, 45)"];   32 [pos="(35.2644, 45)"];   33   [pos="(0,        -30)"];    34 [pos="(-30,     0)"];
35 [pos="(30,       0)"];    36 [pos="(0,        30)"];   37 [pos="(0,       0)"];    0--1 [edge_id=4,      weight=1];
0--2  [edge_id=10, weight=1];     0--3  [edge_id=16, weight=1];     0--4  [edge_id=22, weight=1];     1--2  [edge_id=11, weight=1];     1--3  [edge_id=17, weight=1];
1--5  [edge_id=28, weight=1];     1--6  [edge_id=34, weight=1];     1--7  [edge_id=40, weight=1];     2--4  [edge_id=23, weight=1];     2--5  [edge_id=29, weight=1];
2--8  [edge_id=46, weight=1];     2--10  [edge_id=58, weight=1];    3--4  [edge_id=24, weight=1];     3--6  [edge_id=35, weight=1];     3--9  [edge_id=52, weight=1];
3--11  [edge_id=64, weight=1];    4--8  [edge_id=47, weight=1];     4--9  [edge_id=53, weight=1];     4--12  [edge_id=70, weight=1];    5--7  [edge_id=41, weight=1];
5--10  [edge_id=59, weight=1];    5--14  [edge_id=80, weight=1];    5--16  [edge_id=92, weight=1];    6--7  [edge_id=42, weight=1];     6--11  [edge_id=65, weight=1];
6--15  [edge_id=86, weight=1];    6--17  [edge_id=98, weight=1];    7--13  [edge_id=76, weight=1];    7--14  [edge_id=81, weight=1];    7--15  [edge_id=87, weight=1];
8--10  [edge_id=60, weight=1];    8--12  [edge_id=71, weight=1];    8--20  [edge_id=112, weight=1];   8--22  [edge_id=124, weight=1];   9--11  [edge_id=66, weight=1];
9--12  [edge_id=72, weight=1];    9--21  [edge_id=118, weight=1];   9--23  [edge_id=130, weight=1];   10--16  [edge_id=93, weight=1];   10--18  [edge_id=104, weight=1];
10--20  [edge_id=113, weight=1];  11--17  [edge_id=99, weight=1];   11--19  [edge_id=108, weight=1];  11--21  [edge_id=119, weight=1];  12--22  [edge_id=125, weight=1];
12--23  [edge_id=131, weight=1];  12--24  [edge_id=136, weight=1];  13--14  [edge_id=82, weight=1];   13--15  [edge_id=88, weight=1];   13--25  [edge_id=140, weight=1];
14--16  [edge_id=94, weight=1];   14--25  [edge_id=141, weight=1];  14--28  [edge_id=158, weight=1];  15--17  [edge_id=100, weight=1];  15--25  [edge_id=142, weight=1];
15--29  [edge_id=164, weight=1];  16--18  [edge_id=105, weight=1];  16--26  [edge_id=146, weight=1];  16--28  [edge_id=159, weight=1];  17--19  [edge_id=109, weight=1];
17--27  [edge_id=152, weight=1];  17--29  [edge_id=165, weight=1];  18--20  [edge_id=114, weight=1];  18--26  [edge_id=147, weight=1];  19--21  [edge_id=120, weight=1];
19--27  [edge_id=153, weight=1];  20--22  [edge_id=126, weight=1];  20--26  [edge_id=148, weight=1];  20--31  [edge_id=176, weight=1];  21--23  [edge_id=132, weight=1];
21--27  [edge_id=154, weight=1];  21--32  [edge_id=182, weight=1];  22--24  [edge_id=137, weight=1];  22--30  [edge_id=170, weight=1];  22--31  [edge_id=177, weight=1];
23--24  [edge_id=138, weight=1];  23--30  [edge_id=171, weight=1];  23--32  [edge_id=183, weight=1];  24--30  [edge_id=172, weight=1];  25--28  [edge_id=160, weight=1];
25--29  [edge_id=166, weight=1];  25--33  [edge_id=188, weight=1];  26--28  [edge_id=161, weight=1];  26--31  [edge_id=178, weight=1];  26--34  [edge_id=194, weight=1];
27--29  [edge_id=167, weight=1];  27--32  [edge_id=184, weight=1];  27--35  [edge_id=200, weight=1];  28--33  [edge_id=189, weight=1];  28--34  [edge_id=195, weight=1];
29--33  [edge_id=190, weight=1];  29--35  [edge_id=201, weight=1];  30--31  [edge_id=179, weight=1];  30--32  [edge_id=185, weight=1];  30--36  [edge_id=206, weight=1];
31--34  [edge_id=196, weight=1];  31--36  [edge_id=207, weight=1];  32--35  [edge_id=202, weight=1];  32--36  [edge_id=208, weight=1];  33--34  [edge_id=197, weight=1];
33--35  [edge_id=203, weight=1];  33--37  [edge_id=212, weight=1];  34--36  [edge_id=209, weight=1];  34--37  [edge_id=213, weight=1];  35--36  [edge_id=210, weight=1];
35--37  [edge_id=214, weight=1];  36--37  [edge_id=215, weight=1];
}

O resultado result_ll_pos_0.json:

{
    "0": [1.3505670147100813e-12, 180.0],
    "1": [0.0, -150.0],
    "2": [-30.000000000000004, 180.0],
    "3": [29.999999999999996, 180.0],
    "4": [0.0, 150.0],
    "5": [-35.264389682754654, -135.0],
    "6": [35.264389682754654, -135.0],
    "7": [0.0, -120.0],
    "8": [-35.264389682754654, 135.0],
    "9": [35.264389682754654, 135.0],
    "10": [-59.99999999999999, 180.0],
    "11": [59.99999999999999, 180.0],
    "12": [0.0, 120.0],
    "13": [0.0, -90.0],
    "14": [-29.999999999999996, -90.0],
    "15": [29.999999999999996, -90.0],
    "16": [-59.99999999999999, -90.0],
    "17": [59.99999999999999, -90.0],
    "18": [-89.99999999999866, -0.0006610126423928427],
    "19": [90.0, 0.0],
    "20": [-60.00000000000001, 90.0],
    "21": [60.00000000000001, 90.0],
    "22": [-30.000000000000004, 90.0],
    "23": [29.999999999999996, 90.0],
    "24": [0.0, 90.0],
    "25": [0.0, -60.00000000000001],
    "26": [-60.00000000000001, 0.0],
    "27": [60.00000000000001, 0.0],
    "28": [-35.264389682754654, -45.0],
    "29": [35.264389682754654, -45.0],
    "30": [0.0, 60.00000000000001],
    "31": [-35.264389682754654, 45.0],
    "32": [35.264389682754654, 45.0],
    "33": [0.0, -29.999999999999996],
    "34": [-30.000000000000004, 0.0],
    "35": [29.999999999999996, 0.0],
    "36": [0.0, 29.999999999999996],
    "37": [0.0, 0.0]
}

Eu não poderia ter chegado aqui sem a ajuda de @sehe. Como mencionado nos comentários, eu precisava criar uma estrutura point_difference adequada, e então também notei que não estava alternando entre radianos (boost, matemática C++) e graus (GeographicLib) corretamente. Depois disso, suspeitei que, por algum motivo, o delta sendo passado por kamada_kawai_spring_layout para adjust() as posições precisa ser invertido (ou não entendi a intenção da chamada(!) então adicionei pi/180° ao valor que ele envia.

Embora ainda haja coisas aqui que poderiam ser otimizadas (assim como diversas funções supérfluas), o código funcional é o seguinte.

#include "rapidjson/document.h"
#include "rapidjson/ostreamwrapper.h"
#include "rapidjson/writer.h"
#include <GeographicLib/Constants.hpp>
#include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/Intersect.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/graphviz.hpp>
#include <boost/graph/kamada_kawai_spring_layout.hpp>
#include <boost/graph/random_layout.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_expression.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <fstream>
#include <iostream>
#include <random>
#include <numbers>


static constexpr double r_mux  = std::numbers::pi / 180.0;
static constexpr double d_mux  = 180.0 / std::numbers::pi;

namespace GL = GeographicLib;

template <typename RandomNumberGenerator = std::minstd_rand>
class geo_topology : public boost::convex_topology<2> {
  public:
    using point_type            = typename convex_topology<2>::point_type;
    using rand_t                = boost::uniform_01<RandomNumberGenerator, double>;
    using geo_vec               = boost::numeric::ublas::vector<double>;

  private:
    RandomNumberGenerator gen_;
    mutable rand_t        rand_{gen_};
    GL::Geodesic          geodesic_;

  public:
    geo_topology(geo_topology const&)            = delete;
    geo_topology& operator=(geo_topology const&) = delete;

    explicit geo_topology(GL::Geodesic const& g) : geodesic_(g) {}
    geo_topology(RandomNumberGenerator gen) : gen_(std::move(gen)), geodesic_(GL::Geodesic::WGS84()) {}

    public:
        struct pt_difference {
        private:
            geo_vec v;
        public:
            pt_difference() : v(2) {v[0] = 0.0; v[1] = 0.0;}
            double& operator[](std::size_t i) { return v[i]; }
            const double& operator[](std::size_t i) const { return v[i]; }
                        
            pt_difference(double mag, double heading) : v(2) {
                v[0] = cos(heading) * mag; v[1] = sin(heading) * mag;
            }
            
            double magnitude() const {
                return boost::numeric::ublas::norm_2(v);
            }
            
            double heading() const {
                return atan2(v[1], v[0]);
            }

            friend pt_difference operator+(
                const pt_difference& a, const pt_difference& b) {
                return a.v+b.v;
            }

            friend pt_difference& operator+=(pt_difference& a, const pt_difference& b) {
                a.v = a.v += b.v;
                return a;
            }

            friend pt_difference operator-(const pt_difference& a) {
                pt_difference r;
                return r - a;
            }

            friend pt_difference operator-(const pt_difference& a, const pt_difference& b) {
                pt_difference r;
                r.v = a.v - b.v;
                return r;
            }

            friend pt_difference& operator-=(pt_difference& a, const pt_difference& b) {
                a.v-=b.v;
                return a;
            }

            friend pt_difference operator*(
                const pt_difference& a, const pt_difference& b) {
                return a.v*b.v;
            }

            friend pt_difference operator*(const pt_difference& a, double b) {
                throw std::runtime_error("vector * double");
                return a;
            }

            friend pt_difference operator*(double a, const pt_difference& b){
                double m = b.magnitude() * a;
                double d = b.heading();
                return pt_difference(m, d);
            }

            friend pt_difference operator/(
                const pt_difference& a, const pt_difference& b) {
                return a.v / b.v;
            }

            friend double dot(const pt_difference& a, const pt_difference& b) {
                return boost::numeric::ublas::inner_prod(a.v, b.v);
            }
            
        };

    using point_difference_type = typename geo_topology::pt_difference;

    
    point_type random_point() const {
        throw std::runtime_error("random_point used.");
        point_type p;
        p[0] = fmod(rand_(), std::numbers::pi) - 0.5*std::numbers::pi;  // latitude.
        p[1] = fmod(rand_(), 2.0*std::numbers::pi) - std::numbers::pi; // longitude.
        if (std::isnan(p[0]) or std::isnan(p[1])) {
            std::cerr << "nan in random_point" << std::endl;
        }
        return p;
    }

    point_type bound(point_type a) const {
        throw std::runtime_error("bound used.");
        a[0] = fmod(a[0] + 0.5*std::numbers::pi, std::numbers::pi) - 0.5*std::numbers::pi;   // latitude.
        a[1] = fmod(a[1] + std::numbers::pi, 2.0*std::numbers::pi) - std::numbers::pi; // longitude.
        if (std::isnan(a[0]) or std::isnan(a[1])) {
            std::cerr << "nan in bound" << std::endl;
        }
        return a;
    }

    double distance_from_boundary(point_type /*a*/) const {
        throw std::runtime_error("distance_from_boundary used.");
        return 1000.0; // really don't know if this means anything.
    }

    point_type center() const {
        throw std::runtime_error("center used.");
        point_type result;
        result[0] = std::numbers::pi * 0.5; // latitude.
        result[1] = 0.0;
        return result;
    }

    point_type origin() const {
        throw std::runtime_error("origin used.");
        point_type result;
        result[0] = 0.0; // latitude.
        result[1] = 0.0;
        return result;
    }

    pt_difference extent() const {
        throw std::runtime_error("extent used.");
        return pt_difference(geodesic_.EquatorialRadius(), 2 * std::numbers::pi);
    }

    double distance(point a, point b) const {
        double dist = 0.0;
        double r_fwd_azi, r_rev_azi;
        geodesic_.Inverse(a[0], a[1], b[0], b[1], dist, r_fwd_azi, r_rev_azi);
        return dist;
    }

    point move_position_toward(point a, double fraction, point b) const {
        throw std::runtime_error("move_position_toward used.");
        point  c;
        double dist = 0.0;
        double r_fwd_azi, r_rev_azi;
        geodesic_.Inverse(a[0], a[1], b[0], b[1], dist, r_fwd_azi, r_rev_azi);
        GL::GeodesicLine gl = geodesic_.Line(a[0], a[1], r_fwd_azi);
        gl.Position(fraction * dist, c[0], c[1]);
        return c;
    }

    pt_difference difference(point a, point b) const { // a-b
        double           r_rev_azi;
        double           mag,heading;
        geodesic_.Inverse(a[0], a[1], b[0], b[1], mag, heading, r_rev_azi);
        return pt_difference(mag, r_mux * heading);
    }

    point adjust(point a, pt_difference delta) const {
        point result(a); // point_difference = dist/fwd_azi
        double mag = delta.magnitude();
        double dir = delta.heading() + std::numbers::pi;
        GL::GeodesicLine gl = geodesic_.Line(a[0], a[1], d_mux * dir);
        gl.Position(mag, result[0], result[1]);
        return result;
    }

    point pointwise_min(point a, point b) const {
        throw std::runtime_error("pointwise_min used.");
        point  result;
        double a_dist, b_dist = 0.0;
        double r_fwd_azi, r_rev_azi;
        geodesic_.Inverse(0.5 * std::numbers::pi, 0.0, a[0], a[1], a_dist, r_fwd_azi, r_rev_azi);
        geodesic_.Inverse(0.5 * std::numbers::pi, 0.0, b[0], b[1], b_dist, r_fwd_azi, r_rev_azi);
        result = a_dist <= b_dist ? a : b;
        return result;
    }

    point pointwise_max(point a, point b) const {
        throw std::runtime_error("pointwise_max used.");
        point  result;
        double a_dist, b_dist = 0.0;
        double r_fwd_azi, r_rev_azi;
        geodesic_.Inverse(0.5 * std::numbers::pi, 0.0, a[0], a[1], a_dist, r_fwd_azi, r_rev_azi);
        geodesic_.Inverse(0.5 * std::numbers::pi, 0.0, b[0], b[1], b_dist, r_fwd_azi, r_rev_azi);
        result = a_dist >= b_dist ? a : b;
        return result;
    }

    double norm(pt_difference delta) const {
        return delta.magnitude();
    }

    double volume(pt_difference delta) const {
        throw std::runtime_error("volume used.");
        return delta[1];
    }
};

namespace boost {
    // for graphviz output of node position
    static std::ostream& operator<<(std::ostream& os, geo_topology<>::point const& p) { return os << "(" << p[0] << ", " << p[1] << ")"; }
    static std::istream& operator>>(std::istream& is, geo_topology<>::point&) { is.setstate(std::ios::failbit); return is; } // NOT IMPLEMENTED
} // namespace boost

struct ExampleSpherical {

    static constexpr double radius = 6371008.7714;
    static constexpr double q_circ = radius * std::numbers::pi * 0.5;

    using Topology           = geo_topology<>;
    using Point              = Topology::point_type;
    using EdgeWeightProperty = boost::property<boost::edge_weight_t, double>;

    using Graph = boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS, boost::no_property,
                                        boost::property<boost::edge_index_t, int, EdgeWeightProperty>>;
    using V = boost::graph_traits<Graph>::vertex_descriptor;
    using E = boost::graph_traits<Graph>::edge_descriptor;

    size_t const depth;
    double       line_length;
    using PosVector   = std::vector<Point>;
    using vertex_id_t = unsigned long; // std::size_t;

    std::map<vertex_id_t, std::set<vertex_id_t>> network_;
    std::map<vertex_id_t, Point>                 initial_positions_;
    GL::Geodesic                                 geo;
    Topology                                     topo;

    ExampleSpherical(size_t d) : depth(d), geo(GL::Geodesic::WGS84()), topo(geo) {
        line_length = q_circ / pow(3.0, double(depth + 1));
    }

    void load_positions(std::string const filename) {
        //  Load latitude/longitude positions.
        auto                start = std::chrono::system_clock::now();
        rapidjson::Document jpos;
        std::ifstream       pos_file(filename); // dict of node: pos
        std::string         pos_str(std::istreambuf_iterator<char>(pos_file), {});
        jpos.Parse(pos_str.c_str());
        if (!jpos.HasParseError()) {
            for (auto p_it = jpos.MemberBegin(), n_end = jpos.MemberEnd(); p_it != n_end; ++p_it) {
                auto  node = std::stoul(p_it->name.GetString());
                Point p;
                p[0] = p_it->value[0].GetDouble();
                p[1] = p_it->value[1].GetDouble();
                initial_positions_.emplace(node, std::move(p));
            }
        }
        for (auto& [a_id, a_pos] : initial_positions_) {
            for (auto& [b_id, b_pos] : initial_positions_) {
                if (a_id != b_id) {
                    double mag = topo.difference(a_pos, b_pos).magnitude();
                    if (mag < 0.01) {
                        std::cout << "found overlapping nodes " << a_pos << " " << b_pos << std::endl;
                        exit(0);
                    }
                }
            }
        }
        std::cout << "Optimal Node Distance:" << line_length << std::endl;
        auto end = std::chrono::system_clock::now();
        auto ms  = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
        std::cout << "position load took: " << ms << "ms" << std::endl;
    }

    void load_network(std::string const filename) {
        auto                start = std::chrono::system_clock::now();
        rapidjson::Document jnet;
        std::ifstream net_file(filename); // dict with list "nodes" = each= eg {"id"=0}, and list links=eg
                                          // {"source": 0, "target":2}
        std::string net_str(std::istreambuf_iterator<char>(net_file), {});
        jnet.Parse(net_str.c_str());
        if (!jnet.HasParseError()) {
            for (auto&& n : jnet["nodes"].GetArray())
                network_.insert({n.GetObject()["id"].GetUint64(), {}});

            for (auto&& l : jnet["links"].GetArray()) {
                auto ia = l.GetObject()["source"].GetUint64();
                auto ib = l.GetObject()["target"].GetUint64();
                network_[ia].insert(ib);
                network_[ib].insert(ia);
            }
        }
        auto end = std::chrono::system_clock::now();
        auto ms  = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
        std::cout << "network load took: " << ms << "ms\n";
    }

    void validate_graph() const {
        assert(network_.size() == initial_positions_.size());

        std::set<vertex_id_t> error_nodes;
        for (auto& [id, dd] : network_) {
            if (!initial_positions_.contains(id))
                error_nodes.insert(id); // no position
            if (dd.empty())    // no edges
                error_nodes.insert(id);
            for (auto& q : dd)
                if (!network_.contains(q)) // missing neighbour
                    error_nodes.insert(id);
        }

        if (!error_nodes.empty())
            throw std::runtime_error("files did not reconcile correctly");
    }

    void load_and_validate() {
        load_network("graph_" + std::to_string(depth) + ".json");
        load_positions("graph_ll_pos_" + std::to_string(depth) + ".json");
        validate_graph();
    }

    void save_positions(PosVector& final_pos) const {
        auto start = std::chrono::system_clock::now();
        rapidjson::Document jpos;
        jpos.SetObject();
        auto& alloc = jpos.GetAllocator();
        for (auto& [id, _] : initial_positions_) {
            rapidjson::Value array(rapidjson::kArrayType);
            array.PushBack(final_pos[id][0], alloc);
            array.PushBack(final_pos[id][1], alloc);
            auto j_id = rapidjson::Value(std::to_string(id).c_str(), alloc);
            jpos.AddMember(j_id, array, alloc);
        }
        std::ofstream ofs{"result_ll_pos_" + std::to_string(depth) + ".json"};
        if (!ofs.is_open()) {
            std::cerr << "Could not open file for writing!\n";
        }
        rapidjson::OStreamWrapper                    osw{ofs};
        rapidjson::Writer<rapidjson::OStreamWrapper> writer2{osw};
        jpos.Accept(writer2);

        auto end = std::chrono::system_clock::now();
        auto ms  = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
        std::cout << "position save took: " << ms << "ms\n";
    }

    struct kamada_kawai_done {

        bool operator()(double delta_p, V /*p*/, Graph const& /*g*/, bool global) {
            counter--;
            if (global) {
                auto diff = fabs(last_delta - delta_p);
                std::cout << counter << "; delta_p: " << delta_p << " diff: " << diff << std::endl;
                last_delta = delta_p;
                return fabs(diff) < 0.01 or counter < 1;
            } else {
                if (counter % 100 == 0) {
                    std::cout << counter << "; " << last_delta << std::endl;
                }
                return delta_p < 0.01 or counter < 1;
            }
        }
        double last_delta = 0;
        signed counter = 10000;
    };

    void evaluate(signed cycles) const {
        // std::vector<unsigned long> vertices_(network_.size());
        // for (auto [id, _] : network_)
        //     vertices_.push_back(id);

        Graph     g(network_.size());
        PosVector position(num_vertices(g));
        auto      wmap = get(boost::edge_weight, g);
        auto      eidx = get(boost::edge_index, g);

        for (size_t edge_idx = 0; auto [s, dd] : network_) {
            for (auto d : dd) {
                assert(d != s);
                assert(d < num_vertices(g));
                assert(s < num_vertices(g));
                auto [e, not_duplicate] = add_edge(s, d, g);

                if (not_duplicate) {
                    eidx[e] = edge_idx++; // required to map edges to [0, num_edges(g))
                    wmap[e] = 1.0;        // weight must all be non-negative
                }
            }
            position[s] = initial_positions_.at(s); // throws if not found.
        }
        auto position_map = make_iterator_property_map(position.begin(), get(boost::vertex_index, g));

        {
            // save graph-viz
            boost::dynamic_properties dp;
            dp.property("node_id", get(boost::vertex_index, g));
            dp.property("pos", position_map);
            dp.property("edge_id", eidx);
            dp.property("weight", wmap);

        }

        kamada_kawai_done condition{ line_length, cycles};
        kamada_kawai_spring_layout(g, position_map, wmap, topo, boost::edge_length(line_length), std::ref(condition));
        save_positions(position);
    }
};

int main() {
    ExampleSpherical p(1);
    p.load_and_validate();
    p.evaluate(10000);
}

A rede inicial se parece com isso:

A rede processada/kamada_kawai_spring_layouted é a seguinte:

Esta é uma excelente melhoria. Definitivamente pode haver ajustes, melhorias adicionais e eficiências que poderiam ser implementadas - e eu adoraria ouvir sobre elas - mas a pergunta está respondida. Os Kudos vão para sehe, no entanto.