Konchog's questions

我正在尝试研究如何在 google-sheets 中编写一个直观的公式，给定两个表（布尔值），其维度为AxB，CxB，一个表CxA根据布尔值生成相关表OR。

可以使用什么公式来生成表格CxA？我绝对希望使用单个单元格公式来填充表格。

我一直在考虑使用map、byrow或，但每次开始时我都会迷失在矩阵中bycol。scan

下表重点介绍了以下情况：R对于 来说为真[a,e]； [1,2,3,6]对于 来说为真a；并且[3,4,6]对于 e 来说为真：- 因此[1,2,3,4,6]对于 来说为真R。

我有一个（不可逆）函数 ak([u,v,w])

这将获取单位八面体表面上的一个点（p：使得 |u|+|v|+|w| = 1）并返回单位球体表面上的一个点。该函数并不完美，但目的是保持点之间的距离为等积。

我正在考虑使用 SciPy 最小化来提供数值逆，但我无法理解它。

输入：球面 pt [x,y,z]，
输出八面体 pts [u,v,w]，使得 ak([u,v,w])=[x,y,z]

我的函数 ak() 定义如下：

def ak(p):
    # Convert point on octahedron onto the sphere.
    # Credit to Anders Kaseorg: https://math.stackexchange.com/questions/5016695/
    # input:  oct_pt is a Euclidean point on the surface of a unit octagon.
    # output: UVW on a unit sphere.
    a = 3.227806237143884260376580641604959964752197265625  # 𝛂 - vis. Kaseorg.
    p1 = (np.pi * p) / 2.0
    tp1 = np.tan(p1)
    xu, xv, xw = tp1[0], tp1[1], tp1[2]
    u2, v2, w2 = xu ** 2, xv ** 2, xw ** 2
    y0p = xu * (v2 + w2 + a * w2 * v2) ** 0.25
    y1p = xv * (u2 + w2 + a * u2 * w2) ** 0.25
    y2p = xw * (u2 + v2 + a * u2 * v2) ** 0.25
    pv = np.array([y0p, y1p, y2p])
    return pv / np.linalg.norm(pv, keepdims=True)

此函数基于我在 Math StackExchange 上发表的一篇文章。

有什么提示吗？

我正在尝试让 kamada kawai 弹簧布局在椭圆体表面上工作。我编写了一个 convex_topology 专业化，它实现了 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]}

更新

Sehe 的以下贡献解决了我代码中的许多错误 - 特别是与我对 boost 的 kamada_kawai_spring_layout 的理解不佳有关

拓扑“范数”方法应返回 point_difference 的绝对值（即 {distance, bearing} 对偶）。point_difference 增量不应以负值到达，但它发生在 kamada_kawai 内部的某个地方。

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

我现在看到弹簧值比初始值更能扭曲位置（这在网络的更高分辨率上显示，但对问题没有影响。请注意，其余顶点与原始顶点相比没有变化：

kamada_kawai 看起来应该是这样的（虽然这是通过使用 python 上的 networkx 绘制的，适合单位球体，而不是椭圆体 - 虽然 networkx 很好，但它太慢、占用大量内存，并且不能支持任意拓扑）。

这显示了初始位置..

我正在使用 Python 脚本将输出转换为 KMZ 以实现可视化。

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}')