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:
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).
This shows the initial positions..
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:
#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];
}
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]
}
- How portable is code with #pragma optimize?
- Make warn_unused_result applied to all function with GCC
- Multiple threads and CPU cache
- Linking error vs. compilation error
- how to compile multiple instances of a source file linking to it different files.o
- Background compile in VS2010 and C/C++
- What Turing complete macro languages are available for C/C++
- How to use Cython to compile Python 3 into C
- Are forward declarations needed when the typedef declaration is done?
- How to assign real and imaginary parts of complex variables in C
- Strict aliasing between Complex numbers and real numbers in C
- How to CORRECTLY install gsl library in Linux?
- Using = operator on atomic variable?
- Unknown conversion types in this program I'm meant to compile
- Iterator in C language
- Treating __func__ as a string literal instead of a predefined identifier
- c - strcpy with struct string returning error
- Segmentation fault when executing a Python script in a C program?
- Questions about C Function Prototypes and Compilation
- What is the purpose of "-Wa,option" in GCC?
- Do we have c99 subflags
- netcat gcc compile option so that IDA pro can show function name
- What is the correct way to free output buffer after call to openssl::i2d_X509?
- How can I link C++ files to a C program?
- Why doesn't this program in C compile? Error: undefined reference to `i2c_smbus_read_byte_data'
- integer promotion isn't taken into account
- Why do I get a SEGFAULT when I try to sort an array of integers using K & R code examples?
- Are C preprocessor statements a part of the C language?
- FFMPEG Api conversion from YUV420P to RGB produces strange output
- Why is the tree constructed this way incorrect?