c++graphboostkamada-kawai

Boost kamada_kawai_spring_layout using Topology derived from GeographicLib ... failing


I am trying to get a kamada kawai spring layout to work on the surface of an ellipsoid. I have written a convex_topology specialisation that implements GeographicLib's Geodesic functions, and then I applied a network (stored in json) along with initial positions (likewise stored).

However my code fails. I do not know if it's because the topology is badly defined, or if my implementation of Boost' kamada kawai layout is missing something it needs, or if I am just expecting way too much from everything.

The current topology is storing points as latitude/longitude doubles, and point_difference as a distance / bearing pair, with min/max using distance from the North Pole. It may be better to use the xyz coordinates that Geodesic provides (surface on the ellipsoid, but when I tried that it failed also.

Here is the code I am using..

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

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

Update

Sehe's contribution below resolved many bugs in my code - especially in relation to my poor understanding of boost's kamada_kawai_spring_layout

Topological 'norm' method should return the absolute of the point_difference (which is a {distance, bearing} couplet). The point_difference delta should not be arriving with a negative value, but it's happening somewhere inside kamada_kawai.

    double norm(point_difference delta) const {
        return abs(delta[0]);
    }

I am now seeing the spring values are distorting the positions more than the initial values (this is shown on a higher resolution of the network, but it makes no difference to the issue. Note that the rest of the vertices have not changed from the original: eek

This is what kamada_kawai should be looking somewhat like (though this was drawn fitting a unit sphere, rather than an ellipsoid, via using networkx on python - and, while networkx is good, it's too slow, memory hungry, and cannot support arbitrary topologies). better

This shows the initial positions.. enter image description here

I am using a python script to convert the outputs to KMZ for visualisation..

import json
import simplekml
from utility.spherical import Spherical


def add_line(kz, pts, idx, name=None, col=(154, 154, 154), alpha=0.5, width=4):
    tpt = [pts[i] for i in idx]
    if len(tpt[0]) == 3:
        tpt = [sp.xyz_ll(tf) for tf in tpt]
    lo = [i[1] for i in tpt]
    mim = max(lo) - min(lo)
    if mim > 179.9999:
        path = [(o if o >= 0 else o + 360.0000, a) for (a, o) in tpt]
    else:
        path = [(o, a) for (a, o) in tpt]
    lstr = kz.newlinestring(name=name)
    lstr.coords = path
    lstr.style.linestyle.width = width


if __name__ == '__main__':
    sp = Spherical()
    for depth in [2]:
        pos_f = open(f'result_ll_pos_{depth}.json', 'r')
        pos_n = json.load(pos_f)
        pos_f.close()
        pos = dict()
        for k, v in pos_n.items():
            key = int(k)
            pos[key] = tuple(v)
        net_f = open(f'graph_{depth}.json', 'r')
        net = json.load(net_f)
        net_f.close()
        edges = net['links']
        kml = simplekml.Kml(open=1)
        for edge in edges:
            add_line(kml, pos, [edge['source'], edge['target']])
        kml.save(f'kkp_l_{depth}.kml')
        print(f'kml file generated {depth}')

Solution

  • So after reviewing, I came up with a version that completes.

    There was the bug that inserted 0 as a neighbour to any node (the initializer {{}, {}} creates set with the single element 0, because both {} value-initialze to that value).

    Also, the network stores neighbours in both directions, but since the graph is undirected anyways, this results in all edges being duplicated. If this is not the intention, simply replace vecS with setS for the edges and the duplicates will auto-filter out.

    Then, the edge index property is required, but wasn't populated correctly. They were all set to 1 (1.0, so it looks like that was /intended/ as the weight).

    Neither was the weight. They were all defaulted to 0.

    I think it's even better to not start with last_delta = 0.

    Lastly, it's probably fine, but I wanted to make sure that the last_delta is actually retained, by passing the 'done' condition by reference.

    Finally, I added a write_graphviz step to render the input as parsed for inspection:

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

    This prints e.g.:

    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
    

    The dot graph:

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

    enter image description here

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