* vd, edge PQ vd 3
* l
*
- * By symmetry, this calculation gives the same answer
- * with R and S exchanged. Looking at the projection in
- * the RMS plane:
+ * (The dimensions of this are those of F_vd.)
+ *
+ * By symmetry, this calculation gives the same answer with R and S
+ * exchanged. Looking at the projection in the RMS plane:
*
*
* S'
*
* In practice to avoid division by zero we'll add epsilon to l^3
* and the huge energy ought then to be sufficient for the model to
- * avoid being close to R=S. */
+ * avoid being close to R=S.
+ */
-static double hypotD(const double p[], const double q[]) {
+double hypotD(const double p[D3], const double q[D3]) {
int k;
double pq[D3];
gsl_vector v;
return gsl_blas_snrm2(&v);
}
+double hypotD2(const double p[D3], const double q[D3]) {
+ double d2= 0;
+ K d2= ffsqa(p[k] - q[k], d2);
+ return d2;
+}
+
#ifdef FP_FAST_FMA
# define ffsqa(factor,term) fma((factor),(factor),(term))
#else
-# define ffsqu(factor,term) ((factor)*(factor)+(term))
+# define ffsqa(factor,term) ((factor)*(factor)+(term))
#endif
static const l3_epsison= 1e-6;
K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
K mprime[k]= (vertices[ri][k] + vertices[si][k]) * 0.5;
b= hypotD(vertices[pi], vertices[qi]);
- d2= 0; K d2= ffsqa(m[k] - mprime[k], d2);
+ d2= hypotD2(m, mprime);
l= hypotD(vertices[ri][k] - vertices[si][k]);
l3 = l*l*l + l3_epsilon;
}
-struct VertexIndexMap;
-
-namespace boost {
- struct property_traits<VertexIndexMap> {
- // Concept Readable Property Map:
- typedef int value_type, reference, key_type;
- category
-class Boost
- }};
-
-void single_source_shortest_paths(int v1,
- const double edge_weights[/*f*/],
- ) {
+static void single_source_shortest_paths(int v1,
+ const double edge_weights[/*f*/],
+ double vertex_distances[/*v*/]) {
boost::dijkstra_shortest_paths
(g, v1,
weight_map(edge_weights).
vertex_index_map(identity_property_map()).
-
+ distance_map(vertex_distances));
+}
-void all_pairs_shortest_paths(const Layout *g) {
-
+double graph_layout_energy(const Layout *g) {
+ /* For each (vi,vj) computes shortest path pij = vi..vj along edges,
+ * and actual distance dij = vi-vj.
+ *
+ * Energy contribution is proportional to
+ *
+ * pij - dij
+ * --------- . Sigma length
+ * 2 e member (edges(vi) union edges(vj) e
+ * dij
+ */
+ double edge_weights[N*V6], vertex_distances[N];
+ int v1, e, f;
+
+ FOR_VEDGE_X(v1,e,
+ f= v1 | e << ESHIFT,
+ edge_weights[f]= NaN)
+ edge_weights[f]= hypotD(g.v[v1], g.v[v2]);
+
FOR_VERTEX(v1) {
+ single_source_shortest_paths(v1, edge_weights, vertex_distances);
+ FOR_VERTEX(v2) {
+
+
+
+
+ { ) {
+
0);
extern int edge_end2(unsigned v1, int e);
#define EDGE_END2 edge_end2
+#define FOR_VEDGE_X(v1,e,init,otherwise) \
+ FOR_VPEDGE((v1),(e)) \
+ if (((v2)= EDGE_END2((v1),(e)), \
+ (init), \
+ (v2)) < 0) { otherwise } else
+
+#define NOTHING ((void)0)
+
#define FOR_VEDGE(v1,e) \
- FOR_VPEDGE((v1),(e))
- if (((v2)= EDGE_END2((v1),(e))) < 0) ; else
+ FOR_VEDGE_X(v1,e,NOTHING,NOTHING)
#define FOR_EDGE(v1,e,v2) \
FOR_VERTEX((v1)) \
double v[N][D3];
} Layout;
+double hypotD(const double p[D3], const double q[D3]);
+double hypotD2(const double p[D3], const double q[D3]);
+
#endif /*MGRAPH_H*/