使用从 GeographicLib 派生的拓扑来增强 kamada_kawai_spring_layout ...失败

Question

我正在尝试让 kamada kawai 弹簧布局在椭球体表面上工作。 我编写了一个凸面拓扑专业化，它实现了 GeographicLib 的测地线函数，然后应用了一个网络（存储在 json 中）以及初始位置（同样存储）。

但是我的代码失败了。我不知道是否是因为拓扑定义错误，或者我的 Boost' kamada kawai 布局的实现缺少它需要的东西，或者我只是对一切期望太多。

当前拓扑将点存储为纬度/经度双倍，将 point_difference 存储为距离/方位对，最小/最大使用距北极的距离。 使用Geodesic提供的xyz坐标可能会更好（椭球体上的表面，但是当我尝试时它也失败了。

这是我正在使用的代码..

#include <algorithm>
#include <cassert>
#include <chrono>
#include <cmath>
#include <fstream>
#include <getopt.h>
#include <iostream>
#include <sstream>
#include <streambuf>
#include <string>
#include <unordered_map>
#include <vector>
#include <random>
#include <numbers>

#include <GeographicLib/Intersect.hpp>
#include <GeographicLib/Geodesic.hpp>
#include <GeographicLib/Constants.hpp>

#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/kamada_kawai_spring_layout.hpp>
#include <boost/graph/random_layout.hpp>

#include "rapidjson/document.h"
#include "rapidjson/istreamwrapper.h"
#include "rapidjson/writer.h"
#include "rapidjson/stringbuffer.h"
#include "rapidjson/ostreamwrapper.h"

using namespace std;
using namespace rapidjson;
using namespace GeographicLib;

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

private:
    shared_ptr< RandomNumberGenerator > gen_ptr;
    shared_ptr< rand_t > rand;
    Geodesic wgs84;

public:
    explicit geo_topology(Geodesic& g) : gen_ptr(new RandomNumberGenerator) , rand(new rand_t(*gen_ptr)),wgs84(g) {}
    geo_topology(RandomNumberGenerator& gen): gen_ptr(), rand(new rand_t(gen)),wgs84(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 (isnan(p[0]) or 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 (isnan(a[0]) or 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] = wgs84.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;
        wgs84.Inverse(a[0], a[1], b[0], b[1], dist, r_fwd_azi, r_rev_azi);
        if (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;
        wgs84.Inverse(a[0], a[1], b[0], b[1], dist, r_fwd_azi, r_rev_azi);
        GeodesicLine gl = wgs84.Line(a[0], a[1], r_fwd_azi);
        gl.Position(fraction * dist, c[0], c[1]);
        if (isnan(c[0]) or 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;
        wgs84.Inverse(a[0], a[1], b[0], b[1], result[0], result[1], r_rev_azi);
        if (isnan(result[0]) or 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 (isnan(delta[0]) or 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 {
            GeodesicLine gl = wgs84.Line(a[0], a[1], delta[1]);
            gl.Position(delta[0], result[0], result[1]);
            if (isnan(result[0]) or 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;
        wgs84.Inverse(90.0, 0.0, a[0], a[1], a_dist, r_fwd_azi, r_rev_azi);
        wgs84.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 (isnan(result[0]) or 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;
        wgs84.Inverse(90.0, 0.0, a[0], a[1], a_dist, r_fwd_azi, r_rev_azi);
        wgs84.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 (isnan(result[0]) or isnan(result[1])) {
            std::cerr << "nan in pointwise_max" << std::endl;
        }
        return result;
    }

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

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

};

struct ExampleSpherical {
    
    static constexpr double r_mux = numbers::pi / 180.0;
    static constexpr double d_mux = 180.0 / numbers::pi;
    static constexpr double radius = 6371008.7714;
    static constexpr double q_circ = radius * 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::vecS, boost::vecS, boost::undirectedS,                 //
                                        boost::property<boost::vertex_index_t, int>,                  //
                                        boost::property<boost::edge_index_t, int, EdgeWeightProperty> //
                                        >;
    size_t depth;
    double line_length;
    std::vector<unsigned long> vertices;
    using PosVector = std::vector<Point>;
    using vertex_id_t = unsigned long; //std::size_t;

    map<vertex_id_t,set<vertex_id_t>> net;
    map<vertex_id_t,Point> pos;
    Geodesic geo;
    Topology topo;

    ExampleSpherical(size_t d) : depth(d),geo(Geodesic::WGS84()),topo(geo) {
        line_length = q_circ / pow(3.0, double(depth+1));
    }
    
    void load_positions(const std::string filename) {
    //  Load latitude/longitude positions.
        auto start = std::chrono::system_clock::now();
        Document jpos;
        ifstream pos_file(filename);   // dict of node: pos
        string pos_str((istreambuf_iterator<char>(pos_file)),istreambuf_iterator<char>());
        jpos.Parse(pos_str.c_str());
        if (!jpos.HasParseError()) {
            auto n_end = jpos.MemberEnd();
            for (auto p_it = jpos.MemberBegin(); p_it != n_end; ++p_it) {
                unsigned long node = std::stoul (p_it->name.GetString());
                Point p;
                p[0] = p_it->value[0].GetDouble();
                p[1] = p_it->value[1].GetDouble();
                pos.insert({node,p});
            }
        }
        for (auto& a : pos) {
            for (auto & b : pos) {
                if (a.first != b.first) {
                    auto dx = topo.difference(a.second, b.second);
                    if (dx[0] < line_length) {
                        line_length = dx[0];
                    }
                    if (dx[0] < 0.000001) {
                        cout << "found overlapping nodes";
                        exit(0);
                    }
                }
            }
        }
        cout << "Minimum Node Distance:" << line_length;
        auto end = std::chrono::system_clock::now();
        auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
        cout << "position load took: " << ms << "ms\n";
    }
    void load_network(const std::string filename) {
        auto start = std::chrono::system_clock::now();
        Document jnet;
        ifstream net_file(filename); // dict with list "nodes" = each= eg {"id"=0}, and list links=eg {"source": 0, "target":2}
        string net_str((istreambuf_iterator<char>(net_file)),istreambuf_iterator<char>());
        jnet.Parse(net_str.c_str());
        if (!jnet.HasParseError()) {
            rapidjson::Value::ConstValueIterator n_ite = jnet["nodes"].Begin();
            rapidjson::Value::ConstValueIterator n_end = jnet["nodes"].End();
            do {
                net.insert({n_ite->GetObject()["id"].GetUint64(),{{},{}}});
                n_ite++;
            } while (n_ite != n_end);
            n_ite = jnet["links"].Begin();
            n_end = jnet["links"].End();
            do {
                unsigned long ia = n_ite->GetObject()["source"].GetUint64();
                unsigned long ib = n_ite->GetObject()["target"].GetUint64();
                net[ia].insert(ib);
                net[ib].insert(ia);
                n_ite++;
            } while (n_ite != n_end);
        }
        auto end = std::chrono::system_clock::now();
        auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
        cout << "network load took: " << ms << "ms\n";
    }
    void validate_graph() {
        deque<vertex_id_t> dx;
        for (auto& node : net) {
            if (pos.find(node.first) == pos.end()) {
                dx.push_back(node.first);
            }
        }
        for (auto& node : net) {
            if (node.second.empty()) {
                dx.push_back(node.first);
            }
            for (auto& q : node.second ) {
                if (net.find(q) == net.end()) {
                    dx.push_back(node.first);
                }
            }
        }
        if (!dx.empty()) {
            cerr << "files did not reconcile correctly";
        }
    }
    void load_and_validate() {
        ostringstream graph_file_str;
        ostringstream graph_pos_str;
        graph_file_str << "graph_" << depth << ".json";
        graph_pos_str << "graph_ll_pos_"  << depth << ".json";
        load_network(graph_file_str.str());
        load_positions(graph_pos_str.str());
        validate_graph();
    }
    void save_positions(PosVector& final_pos) {
        ostringstream result_pos_str;
        result_pos_str << "result_ll_pos_"  << depth << ".json";
        auto start = std::chrono::system_clock::now();
        Document jpos;
        jpos.SetObject();
        auto& alloc = jpos.GetAllocator();
        for (auto& node: pos) {
            ostringstream os;
            os << node.first;
            string id = os.str();
            Value array(kArrayType);
            array.PushBack(final_pos[node.first][0],alloc);
            array.PushBack(final_pos[node.first][1],alloc);
            auto j_id = Value(id.c_str(), alloc);
            jpos.AddMember(j_id, array, alloc);
        }
        std::ofstream ofs { result_pos_str.str() };
        if ( !ofs.is_open() ) {
            std::cerr << "Could not open file for writing!\n";
        }
        OStreamWrapper osw { ofs };
        Writer<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();
        cout << "position save took: " << ms << "ms\n";
    }

    struct kamada_kawai_done {
        kamada_kawai_done() : last_delta() {}
        template<typename Graph>
        bool operator()(double delta_p,
                        typename boost::graph_traits<Graph>::vertex_descriptor /*p*/,
                        const Graph& /*g*/,
                        bool global)
        {
            if (global) {
                double diff = last_delta - delta_p;
                if (diff < 0) diff = -diff;
                std::cout << "delta_p: " << delta_p << std::endl;
                last_delta = delta_p;
                return diff < 0.01;
            } else {
                return delta_p < 0.01;
            }
        }
        double last_delta;
    };
    
    void evaluate() {
        load_and_validate();
        vertices.resize(net.size());
        size_t i=0;
        for (auto n: net) {
            vertices[i] = n.first;
            i++;
        }
        Graph g(this->vertices.size());
        PosVector position(num_vertices(g));

        for (auto n: net) {
            auto& s = n.first;
            for (auto d: n.second) {
                if (d != s) {
                    add_edge(s, d, 1.0, g);
                } else {
                    cout << "Weird. source and destinations " << s << " and " << d << " are the same. Skipping that." << endl;
                }
            }
            position[s] = pos[s];
        }
        auto position_map = boost::make_iterator_property_map(position.begin(), get(boost::vertex_index, g));
        boost::kamada_kawai_spring_layout(g, position_map, get(boost::edge_weight, g), topo, boost::edge_length(10.0), kamada_kawai_done());
        save_positions(position);
    }
};

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

graph_0.json

{"directed": false, "multigraph": false, "graph": {}, "nodes": [{"id": 0}, {"id": 1}, {"id": 2}, {"id": 3}, {"id": 4}, {"id": 5}, {"id": 6}, {"id": 7}, {"id": 8}, {"id": 9}, {"id": 10}, {"id": 11}, {"id": 12}, {"id": 13}, {"id": 14}, {"id": 15}, {"id": 16}, {"id": 17}, {"id": 18}, {"id": 19}, {"id": 20}, {"id": 21}, {"id": 22}, {"id": 23}, {"id": 24}, {"id": 25}, {"id": 26}, {"id": 27}, {"id": 28}, {"id": 29}, {"id": 30}, {"id": 31}, {"id": 32}, {"id": 33}, {"id": 34}, {"id": 35}, {"id": 36}, {"id": 37}], "links": [{"source": 0, "target": 2}, {"source": 0, "target": 1}, {"source": 0, "target": 4}, {"source": 0, "target": 3}, {"source": 1, "target": 7}, {"source": 1, "target": 3}, {"source": 1, "target": 5}, {"source": 1, "target": 2}, {"source": 1, "target": 6}, {"source": 2, "target": 4}, {"source": 2, "target": 8}, {"source": 2, "target": 5}, {"source": 2, "target": 10}, {"source": 3, "target": 9}, {"source": 3, "target": 11}, {"source": 3, "target": 6}, {"source": 3, "target": 4}, {"source": 4, "target": 9}, {"source": 4, "target": 8}, {"source": 4, "target": 12}, {"source": 5, "target": 7}, {"source": 5, "target": 14}, {"source": 5, "target": 16}, {"source": 5, "target": 10}, {"source": 6, "target": 15}, {"source": 6, "target": 17}, {"source": 6, "target": 11}, {"source": 6, "target": 7}, {"source": 7, "target": 15}, {"source": 7, "target": 13}, {"source": 7, "target": 14}, {"source": 8, "target": 10}, {"source": 8, "target": 20}, {"source": 8, "target": 12}, {"source": 8, "target": 22}, {"source": 9, "target": 21}, {"source": 9, "target": 12}, {"source": 9, "target": 11}, {"source": 9, "target": 23}, {"source": 10, "target": 20}, {"source": 10, "target": 16}, {"source": 10, "target": 18}, {"source": 11, "target": 17}, {"source": 11, "target": 19}, {"source": 11, "target": 21}, {"source": 12, "target": 22}, {"source": 12, "target": 24}, {"source": 12, "target": 23}, {"source": 13, "target": 14}, {"source": 13, "target": 25}, {"source": 13, "target": 15}, {"source": 14, "target": 28}, {"source": 14, "target": 16}, {"source": 14, "target": 25}, {"source": 15, "target": 25}, {"source": 15, "target": 29}, {"source": 15, "target": 17}, {"source": 16, "target": 26}, {"source": 16, "target": 28}, {"source": 16, "target": 18}, {"source": 17, "target": 27}, {"source": 17, "target": 29}, {"source": 17, "target": 19}, {"source": 18, "target": 26}, {"source": 18, "target": 20}, {"source": 19, "target": 27}, {"source": 19, "target": 21}, {"source": 20, "target": 26}, {"source": 20, "target": 22}, {"source": 20, "target": 31}, {"source": 21, "target": 32}, {"source": 21, "target": 23}, {"source": 21, "target": 27}, {"source": 22, "target": 30}, {"source": 22, "target": 31}, {"source": 22, "target": 24}, {"source": 23, "target": 24}, {"source": 23, "target": 30}, {"source": 23, "target": 32}, {"source": 24, "target": 30}, {"source": 25, "target": 28}, {"source": 25, "target": 29}, {"source": 25, "target": 33}, {"source": 26, "target": 34}, {"source": 26, "target": 31}, {"source": 26, "target": 28}, {"source": 27, "target": 29}, {"source": 27, "target": 35}, {"source": 27, "target": 32}, {"source": 28, "target": 33}, {"source": 28, "target": 34}, {"source": 29, "target": 33}, {"source": 29, "target": 35}, {"source": 30, "target": 32}, {"source": 30, "target": 36}, {"source": 30, "target": 31}, {"source": 31, "target": 36}, {"source": 31, "target": 34}, {"source": 32, "target": 35}, {"source": 32, "target": 36}, {"source": 33, "target": 37}, {"source": 33, "target": 35}, {"source": 33, "target": 34}, {"source": 34, "target": 37}, {"source": 34, "target": 36}, {"source": 35, "target": 37}, {"source": 35, "target": 36}, {"source": 36, "target": 37}]}

和 graph_ll_0.json

{"0": [0.0, 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.00000000000001], "8": [-35.264389682754654, 135.0], "9": [35.264389682754654, 135.0], "10": [-59.99999999999999, 180.0], "11": [59.99999999999999, 180.0], "12": [0.0, 120.00000000000001], "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": [-90.0, 0.0], "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]}

如果我没有花一周时间试图弄清楚这一切，我就不会问这个问题。

Answer 1

经过审阅，我想出了一个完整的版本。

存在将

0

作为邻居插入到任何节点的错误（初始化器

{{}, {}}

使用单个元素

0

创建集合，因为

{}

值都初始化为该值）。

此外，网络还存储两个方向的邻居，但由于图形无论如何都是

undirected

，这会导致所有边都被重复。如果这不是本意，只需将边缘的

vecS

替换为

setS

，重复项就会自动过滤掉。

然后，边索引属性是必需的，但未正确填充。它们都设置为

1

（1.0，所以看起来/打算/作为权重）。

重量也不是。它们都默认为 0。

我认为最好不要从

last_delta = 0

开始。

最后，这可能没问题，但我想通过引用传递“完成”条件来确保

last_delta

实际上被保留。

最后，我添加了一个

write_graphviz

步骤来渲染输入以供检查解析：

住在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;
    std::minstd_rand                             rng{42}; // for reproducibility
    Topology                                     topo{rng};

    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  = add_edge(s, d, g).first;

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

注意

删除随机种子 (42) 以再次获得随机行为！

对于固定种子，打印：

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

点图：

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];
}

结果

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]
}