ʻOhana
Population structure, admixture history, and selection using learning methods.
jade.svg_tree.hpp
1 /* -------------------------------------------------------------------------
2  Ohana
3  Copyright (c) 2015-2020 Jade Cheng (\___/)
4  Jade Cheng <info@jade-cheng.com> (='.'=)
5  ------------------------------------------------------------------------- */
6 
7 #ifndef JADE_SVG_HPP__
8 #define JADE_SVG_HPP__
9 
10 #include "jade.newick.hpp"
11 #include "jade.simplex.hpp"
12 #include "jade.vec2.hpp"
13 
14 namespace jade
15 {
16  ///
17  /// A template for a class that renders newick trees as SVG.
18  ///
19  template <typename TValue>
20  class basic_svg_tree
21  {
22  public:
23  /// The value type.
24  typedef TValue value_type;
25 
26  /// The Newick node type.
27  typedef basic_newick_node<value_type> node_type;
28 
29  /// The vector type.
30  typedef basic_vec2<value_type> vec2_type;
31 
32  ///
33  /// Initializes a new instance of the class based on the specified
34  /// node.
35  ///
36  explicit basic_svg_tree(
37  const node_type & node) ///< The newick node.
38  : _root (node.get_id())
39  , _table ()
40  {
41  _init_table(node);
42  _update_table(node);
43  }
44 
45  ///
46  /// Optimizes the positions of the vertices using the Nelder-Mead
47  /// Simplex method.
48  ///
49  void optimize_positions()
50  {
51  typedef basic_simplex<value_type> simplex_type;
52  typedef typename simplex_type::container_type container_type;
53  typedef typename simplex_type::execute_args execute_args_type;
54  typedef typename simplex_type::options options_type;
55 
56  //
57  // Updates the radians and positions of each node based on a
58  // specified Nelder-Mead parameter set.
59  //
60  const auto import_params = [this](
61  const container_type & params) -> void
62  {
63  auto dst = _table.begin();
64  for (auto src = params.begin(); src != params.end(); ++src)
65  (dst++)->second.radians = *src;
66  _update_table(*_table.find(_root)->second.node);
67  };
68 
69  //
70  // Computes the objective function for the Nelder-Mead method. The
71  // value returned is the sum of the inverses of the distances
72  // between all pairs of distinct nodes.
73  //
74  const auto objfunc = [this, import_params](
75  const container_type & params) -> value_type
76  {
77  import_params(params);
78  auto sum = value_type(0);
79  for (const auto & a0 : _table)
80  for (const auto & b0 : _table)
81  if (a0.first != b0.first)
82  {
83  const auto a = a0.second.position;
84  const auto b = b0.second.position;
85  const auto d = vec2_type::distance_squared(a, b);
86  if (d < 1e-6) // epsilon
87  return std::numeric_limits<value_type>::max();
88  sum += value_type(1) / d;
89  }
90  return sum;
91  };
92 
93  //
94  // Converts a specified number of degrees to radians.
95  //
96  const auto to_radians = [](const value_type deg) -> value_type
97  {
98  static const auto tau = value_type(2.0 * std::acos(-1.0));
99  return deg * tau / value_type(360.0);
100  };
101 
102  const auto n = _table.size();
103 
104  //
105  // Construct the simplex using a unit size of one degree and the
106  // current positions of the nodes.
107  //
108  options_type opts (n);
109  opts.unit = to_radians(1);
110  opts.vertex.clear();
111  for (const auto & p : _table)
112  opts.vertex.push_back(p.second.radians);
113  simplex_type simplex (objfunc, opts);
114 
115  //
116  // Execute the Nelder-Mead method until the simplex has collapsed
117  // or a maximum timeout occurs; execute at least a minimum number
118  // of iterations.
119  //
120  execute_args_type exe_args_1;
121  exe_args_1.max_iterations = 100;
122  exe_args_1.max_seconds = 0.5;
123  simplex.execute(objfunc, exe_args_1);
124 
125  execute_args_type exe_args_2;
126  exe_args_2.min_length = to_radians(0.01) * value_type(n);
127  exe_args_2.max_seconds = 0.5;
128  simplex.execute(objfunc, exe_args_2);
129 
130  //
131  // Import the final parameters from the simplex.
132  //
133  import_params(simplex.get_vertex());
134  }
135 
136  ///
137  /// Returns the SVG scene as a string.
138  ///
139  /// \return The SVG scene as a string.
140  ///
141  std::string str() const
142  {
143  std::ostringstream out;
144  write(out);
145  return out.str();
146  }
147 
148  ///
149  /// Writes the SVG scene to the specified output stream.
150  ///
151  void write(
152  std::ostream & out) ///< The output stream.
153  const
154  {
155  const auto & root = *_table.find(_root)->second.node;
156  metrics_data metrics;
157  _create_metrics(metrics);
158 
159  _write_svg_header(out, metrics);
160  _write_edges(out, metrics, root);
161  _write_nodes(out, metrics, root);
162  _write_svg_footer(out);
163  }
164 
165  ///
166  /// Writes the SVG scene to the specified output file.
167  ///
168  void write(
169  char const * const path) ///< The path to the output file.
170  const
171  {
172  assert(path != nullptr);
173  std::ofstream out (path);
174  if (!out.good())
175  throw error() << "error opening '" << path << "' for writing";
176  write(out);
177  }
178 
179  ///
180  /// Writes the SVG scene to the specified output file.
181  ///
182  void write(
183  const std::string & path) ///< The path to the output file.
184  const
185  {
186  write(path.c_str());
187  }
188 
189  private:
190  // --------------------------------------------------------------------
191  struct node_ex
192  {
193  const node_type * node;
194  value_type length;
195  value_type radians;
196  vec2_type position;
197 
198  // ----------------------------------------------------------------
199  inline node_ex()
200  : node (nullptr)
201  , length ()
202  , radians ()
203  , position ()
204  {
205  }
206  };
207 
208  // --------------------------------------------------------------------
209  struct metrics_data
210  {
211  value_type font_size;
212  value_type large_radius;
213  value_type small_radius;
214  value_type stroke_width;
215  value_type text_offset;
216  value_type scale;
217  value_type view_box_height;
218  value_type view_box_left;
219  value_type view_box_top;
220  value_type view_box_width;
221  };
222 
223  // --------------------------------------------------------------------
224  struct theme
225  {
226  static constexpr double opacity = 0.75;
227  static constexpr char const * circle_fill_color = "black";
228  static constexpr char const * circle_stroke_color = "white";
229  static constexpr char const * text_color = "white";
230  static constexpr char const * text_outline_color = "black";
231  static constexpr char const * line_color = "black";
232  static constexpr char const * font_family = "Times,serif";
233  };
234 
235  typedef std::map<int, node_ex> table_type;
236 
237  // --------------------------------------------------------------------
238  void _create_metrics(metrics_data & metrics) const
239  {
240  typedef std::numeric_limits<value_type> numeric_limits_type;
241 
242  memset(&metrics, 0, sizeof(metrics));
243 
244  auto iter = _table.begin();
245  if (iter == _table.end())
246  return;
247 
248  metrics.scale = value_type(1);
249 
250  const auto is_measured = [](const node_ex & n) -> bool
251  {
252  return n.node->has_length() && n.node->has_name();
253  };
254 
255  auto shortest_edge = numeric_limits_type::quiet_NaN();
256  if (is_measured(iter->second))
257  shortest_edge = iter->second.length;
258 
259  auto min = iter->second.position;
260  auto max = iter->second.position;
261 
262  if (++iter == _table.end())
263  {
264  metrics.large_radius = value_type(100);
265  metrics.scale = value_type(1);
266  }
267  else
268  {
269  while (iter != _table.end())
270  {
271  if (is_measured(iter->second))
272  {
273  const auto shortest_edge_i =
274  iter->second.length;
275 
276  if (shortest_edge_i > value_type(0))
277  if (std::isnan(shortest_edge)
278  || shortest_edge_i < shortest_edge)
279  shortest_edge = shortest_edge_i;
280  }
281 
282  min = vec2_type::min(min, iter->second.position);
283  max = vec2_type::max(max, iter->second.position);
284  ++iter;
285  }
286 
287  if (std::isnan(shortest_edge))
288  {
289  metrics.large_radius = value_type(100);
290  metrics.scale = value_type(1);
291  }
292  else
293  {
294  auto diameter = vec2_type::distance(min, max);
295 
296  if (diameter < value_type(0.000001))
297  {
298  metrics.large_radius = value_type(100);
299  metrics.scale = value_type(1);
300  }
301  else
302  {
303  while (diameter < value_type(100))
304  {
305  metrics.scale *= value_type(10);
306  diameter *= value_type(10);
307  }
308 
309  while (diameter > value_type(1000))
310  {
311  metrics.scale /= value_type(10);
312  diameter /= value_type(10);
313  }
314 
315  metrics.large_radius = metrics.scale *
316  shortest_edge / value_type(3);
317  }
318  }
319  }
320 
321  const auto diameter = metrics.large_radius + metrics.large_radius;
322  min = metrics.scale * min - diameter;
323  max = metrics.scale * max + diameter;
324 
325  metrics.font_size = metrics.large_radius * value_type(0.90);
326  metrics.small_radius = metrics.large_radius * value_type(0.30);
327  metrics.stroke_width = metrics.large_radius * value_type(0.15);
328  metrics.text_offset = metrics.large_radius * value_type(0.30);
329  metrics.view_box_left = min.x;
330  metrics.view_box_top = min.y;
331  metrics.view_box_width = (max - min).x;
332  metrics.view_box_height = (max - min).y;
333  return;
334  }
335 
336  // --------------------------------------------------------------------
337  void _find_edge_coords(
338  const metrics_data & metrics,
339  const node_type & n1,
340  const node_type & n2,
341  vec2_type & p1,
342  vec2_type & p2)
343  const
344  {
345  const auto & n1ex = _table.find(n1.get_id())->second;
346  const auto & n2ex = _table.find(n2.get_id())->second;
347 
348  p1 = metrics.scale * n1ex.position;
349  p2 = metrics.scale * n2ex.position;
350 
351  const auto n = vec2_type::normalize(p2 - p1);
352 
353  p1 += n * (n1ex.node->get_name().empty()
354  ? metrics.small_radius
355  : metrics.large_radius);
356 
357  p2 -= n * (n2ex.node->get_name().empty()
358  ? metrics.small_radius
359  : metrics.large_radius);
360  }
361 
362  // --------------------------------------------------------------------
363  void _init_table(const node_type & node)
364  {
365  static const auto tau = value_type(2.0 * std::acos(-1.0));
366  _init_table(node, 0, tau);
367  }
368 
369  // --------------------------------------------------------------------
370  void _init_table(
371  const node_type & node,
372  const value_type radians,
373  const value_type range)
374  {
375  node_ex ex;
376  ex.node = &node;
377  ex.radians = radians;
378  ex.length = std::max(value_type(0), node.get_length());
379 
380  const auto id = node.get_id();
381  _table[id] = ex;
382 
383  const auto & children = node.get_children();
384  const auto n = children.size();
385  const auto nval = value_type(n);
386  const auto half = value_type(0.5);
387  for (size_t i = 0; i < n; i++)
388  {
389  const auto ival = value_type(i);
390  const auto percent = (half + ival) / nval;
391  const auto radians2 = (range * percent) - (half * range);
392  const auto range2 = range / nval;
393  _init_table(*children[i], radians2, range2);
394  }
395  }
396 
397  // --------------------------------------------------------------------
398  void _update_table(const node_type & node)
399  {
400  _update_table(node, 0, vec2_type());
401  }
402 
403  // --------------------------------------------------------------------
404  void _update_table(
405  const node_type & node,
406  const value_type radians0,
407  const vec2_type & p0)
408  {
409  const auto id = node.get_id();
410  auto & ex = _table[id];
411 
412  const auto rad = radians0 + ex.radians;
413 
414  ex.position = p0 + ex.length * vec2_type(
415  std::cos(rad),
416  std::sin(rad));
417 
418  for (const auto ptr : node.get_children())
419  _update_table(*ptr, rad, ex.position);
420  }
421 
422  // --------------------------------------------------------------------
423  void _write_edges(
424  std::ostream & out,
425  const metrics_data & metrics,
426  const node_type & node)
427  const
428  {
429  const auto stroke_width = metrics.stroke_width;
430 
431  for (const auto child : node.get_children())
432  {
433  vec2_type p1, p2;
434  _find_edge_coords(metrics, node, *child, p1, p2);
435  _write_svg_line(out, p1.x, p1.y, p2.x, p2.y, stroke_width);
436  }
437 
438  for (const auto child : node.get_children())
439  _write_edges(out, metrics, *child);
440  }
441 
442  // --------------------------------------------------------------------
443  void _write_nodes(
444  std::ostream & out,
445  const metrics_data & metrics,
446  const node_type & node)
447  const
448  {
449  const auto stroke_width = value_type(0.5) * metrics.stroke_width;
450  const auto parent_id = node.get_id();
451  const auto & parent_ex = _table.find(parent_id)->second;
452  const auto position = metrics.scale * parent_ex.position;
453 
454  const auto radius = node.get_name().empty() ?
455  metrics.small_radius :
456  metrics.large_radius;
457 
458  out << " <g>\n";
459 
460  _write_svg_circle(out, position.x, position.y,
461  radius, stroke_width);
462 
463  if (!node.get_name().empty())
464  {
465  const auto x = position.x;
466  const auto y = position.y + metrics.text_offset;
467  const auto font_size = metrics.font_size;
468  const auto text = node.get_name().c_str();
469  _write_svg_text(out, x, y, font_size, text, true);
470  _write_svg_text(out, x, y, font_size, text, false);
471  }
472 
473  out << " </g>\n";
474 
475  for (const auto child : node.get_children())
476  _write_nodes(out, metrics, *child);
477  }
478 
479  // --------------------------------------------------------------------
480  static void _write_svg_circle(
481  std::ostream & out,
482  const value_type cx,
483  const value_type cy,
484  const value_type r,
485  const value_type stroke_width)
486  {
487  out << " <circle "
488  << "cx=\"" << cx << "px\" "
489  << "cy=\"" << cy << "px\" "
490  << "r=\"" << r << "px\" "
491  << "style=\""
492  << "fill:" << theme::circle_fill_color << ";"
493  << "fill-opacity:" << theme::opacity << ";"
494  << "stroke:" << theme::circle_stroke_color << ";"
495  << "stroke-width:" << stroke_width << "px\">"
496  << "</circle>\n";
497  }
498 
499  // --------------------------------------------------------------------
500  static void _write_svg_footer(
501  std::ostream & out)
502  {
503  out << "</svg>\n";
504  }
505 
506  // --------------------------------------------------------------------
507  static void _write_svg_header(
508  std::ostream & out,
509  const metrics_data & metrics)
510  {
511  out << "<svg "
512  << "xmlns=\"http://www.w3.org/2000/svg\" "
513  << "version=\"1.1\" "
514  << "viewBox=\""
515  << metrics.view_box_left << " "
516  << metrics.view_box_top << " "
517  << metrics.view_box_width << " "
518  << metrics.view_box_height << "\">\n";
519  }
520 
521  // --------------------------------------------------------------------
522  static void _write_svg_line(
523  std::ostream & out,
524  const value_type x1,
525  const value_type y1,
526  const value_type x2,
527  const value_type y2,
528  const value_type stroke_width)
529  {
530  out << " <line "
531  << "x1=\"" << x1 << "px\" "
532  << "y1=\"" << y1 << "px\" "
533  << "x2=\"" << x2 << "px\" "
534  << "y2=\"" << y2 << "px\" "
535  << "style=\""
536  << "stroke-opacity:" << theme::opacity << ";"
537  << "stroke-linecap:butt;"
538  << "stroke:" << theme::line_color << ";"
539  << "stroke-width:" << stroke_width << "px\">"
540  << "</line>\n";
541  }
542 
543  // --------------------------------------------------------------------
544  static void _write_svg_text(
545  std::ostream & out,
546  const value_type x,
547  const value_type y,
548  const value_type font_size,
549  char const * const text,
550  bool is_outlined)
551  {
552  out << " <text "
553  << "x=\"" << x << "px\" "
554  << "y=\"" << y << "px\" "
555  << "style=\"";
556 
557  if (is_outlined)
558  out << "fill:" << theme::text_outline_color << ";"
559  << "fill-opacity:" << theme::opacity << ";"
560  << "stroke-opacity:" << theme::opacity << ";"
561  << "stroke:" << theme::text_outline_color << ";"
562  << "stroke-width:" << font_size / value_type(7) << "px;";
563  else
564  out << "fill:" << theme::text_color << ";";
565 
566  out << "text-anchor:middle;"
567  << "font-size:" << font_size << "px;"
568  << "font-family:" << theme::font_family << "\">";
569 
570  for (char const * ptr = text; *ptr != '\0'; ptr++)
571  if (std::isprint(*ptr))
572  out << *ptr;
573 
574  out << "</text>\n";
575  }
576 
577  int _root;
578  table_type _table;
579  };
580 }
581 
582 #endif // JADE_SVG_HPP__
jade::basic_svg_tree::value_type
TValue value_type
The value type.
Definition: jade.svg_tree.hpp:24
jade::basic_vec2< value_type >
jade::basic_vec2< value_type >::normalize
static basic_vec2 normalize(const basic_vec2 &v)
Definition: jade.vec2.hpp:137
jade::basic_svg_tree::write
void write(std::ostream &out) const
Writes the SVG scene to the specified output stream.
Definition: jade.svg_tree.hpp:151
jade::basic_newick_node
A template class representing a node from a Newick tree.
Definition: jade.newick.hpp:19
jade::basic_simplex
A template for a class that implements the Nelder-Mead Simplex Method, which attempts to minimize an ...
Definition: jade.simplex.hpp:21
jade::basic_vec2< value_type >::max
static basic_vec2 max(const basic_vec2 &a, const basic_vec2 &b)
Definition: jade.vec2.hpp:117
jade::basic_vec2< value_type >::distance
static value_type distance(const basic_vec2 &a, const basic_vec2 &b)
Definition: jade.vec2.hpp:59
jade::basic_svg_tree::str
std::string str() const
Returns the SVG scene as a string.
Definition: jade.svg_tree.hpp:141
jade::basic_svg_tree::vec2_type
basic_vec2< value_type > vec2_type
The vector type.
Definition: jade.svg_tree.hpp:30
jade::basic_svg_tree::write
void write(char const *const path) const
Writes the SVG scene to the specified output file.
Definition: jade.svg_tree.hpp:168
jade::basic_vec2< value_type >::min
static basic_vec2 min(const basic_vec2 &a, const basic_vec2 &b)
Definition: jade.vec2.hpp:127
jade::basic_vec2< value_type >::distance_squared
static value_type distance_squared(const basic_vec2 &a, const basic_vec2 &b)
Definition: jade.vec2.hpp:69
jade::basic_svg_tree::write
void write(const std::string &path) const
Writes the SVG scene to the specified output file.
Definition: jade.svg_tree.hpp:182
jade::basic_svg_tree
A template for a class that renders newick trees as SVG.
Definition: jade.svg_tree.hpp:21
jade::basic_svg_tree::basic_svg_tree
basic_svg_tree(const node_type &node)
Initializes a new instance of the class based on the specified node.
Definition: jade.svg_tree.hpp:36
jade::basic_error
A template for a class representing an exception thrown from this namespace.
Definition: jade.error.hpp:20
jade::basic_svg_tree::optimize_positions
void optimize_positions()
Optimizes the positions of the vertices using the Nelder-Mead Simplex method.
Definition: jade.svg_tree.hpp:49
jade::basic_svg_tree::node_type
basic_newick_node< value_type > node_type
The Newick node type.
Definition: jade.svg_tree.hpp:27
Copyright © 2015-2020 Jade Cheng