text: //
  // Created by Alan Freitas on 2020-07-05.
  //
  
  #include <algorithm>
  #include <cmath>
  #include <iomanip>
  #include <memory>
  #include <sstream>
  
  #include <matplot/core/axes_object.h>
  #include <matplot/core/axes_type.h>
  #include <matplot/core/figure_type.h>
  
  #include <matplot/util/colors.h>
  #include <matplot/util/common.h>
  #include <matplot/util/concepts.h>
  #include <matplot/util/contourc.h>
  #include <matplot/util/geodata.h>
  
  #include <matplot/axes_objects/bars.h>
  #include <matplot/axes_objects/box_chart.h>
  #include <matplot/axes_objects/circles.h>
  #include <matplot/axes_objects/contours.h>
  #include <matplot/axes_objects/error_bar.h>
  #include <matplot/axes_objects/filled_area.h>
  #include <matplot/axes_objects/function_line.h>
  #include <matplot/axes_objects/labels.h>
  #include <matplot/axes_objects/line.h>
  #include <matplot/axes_objects/matrix.h>
  #include <matplot/axes_objects/network.h>
  #include <matplot/axes_objects/parallel_lines.h>
  #include <matplot/axes_objects/stair.h>
  #include <matplot/axes_objects/string_function.h>
  #include <matplot/axes_objects/surface.h>
  #include <matplot/axes_objects/vectors.h>
  
  #include <matplot/freestanding/histcounts.h>
  
  namespace matplot {
  
      void axes_type::draw() { parent_->draw(); }
  
      void axes_type::touch() { parent_->touch(); }
  
      void axes_type::emplace_object(axes_object_handle obj) {
          if (next_plot_replace_) {
              children_.clear();
              colororder_index_ = 0;
          }
          children_.push_back(obj);
          touch();
      }
  
      void axes_type::run_title_command() {
          include_comment("Axes title");
          if (title_visible_ && !title_.empty()) {
              std::string cmd;
              cmd += "set title \"";
              if (title_enhanced_ && iequals(title_font_weight_, "bold")) {
                  cmd += "{/:Bold ";
              }
              if (title_enhanced_ && iequals(title_font_weight_, "italic")) {
                  cmd += "{/:Italic ";
              }
              cmd += escape(title());
              if (title_enhanced_ && iequals(title_font_weight_, "bold")) {
                  cmd += "}";
              }
              if (title_enhanced_ && iequals(title_font_weight_, "italic")) {
                  cmd += "}";
              }
              if (parent_->backend_->supports_fonts()) {
                  cmd += "\" font \"" + escape(font()) + "," +
                         std::to_string(unsigned(font_size() *
                                                 title_font_size_multiplier_)) +
                         "\"";
                  cmd += " textcolor '" + to_string(title_color_) + "'";
              }
              if (title_enhanced_) {
                  cmd += " enhanced";
              } else {
                  cmd += " noenhanced";
              }
              run_command(cmd);
          } else {
              run_command("unset title");
          }
      }
  
      void axes_type::run_colormap_command() {
          bool any_obj_needs_colormap = false;
          for (const auto &child : children_) {
              if (child->requires_colormap()) {
                  any_obj_needs_colormap = true;
                  break;
              }
          }
          if (any_obj_needs_colormap && !colormap_.empty() &&
              !children_.empty()) {
              std::stringstream ss;
              // limit the number of colors in the palette
              // this is useful for contour plots
              if (max_colors_) {
                  ss << "    set palette maxcolors " << max_colors_ << "\n";
              }
              ss << "    set palette model RGB\n";
              ss << "    set palette defined (\\\n";
              if (max_colors_ == 0) {
                  // put all colormap colors in the palette
                  if (colormap_.size() == 1) {
                      colormap_.emplace_back(colormap_.front());
                  }
                  for (size_t i = 0; i < colormap_.size(); ++i) {
                      if (colormap_[i].size() == 4) {
                          std::array<double, 4> c = to_array<4>(colormap_[i]);
                          std::array<float, 4> c2;
                          std::copy(c.begin(), c.end(), c2.begin());
                          ss << "    " << i << "   \"" << to_string(c2) << "\"";
                      } else {
                          ss << "    " << i << "   " << colormap_[i][0] << " "
                             << colormap_[i][1] << " " << colormap_[i][2];
                      }
                      if (i != colormap_.size() - 1) {
                          ss << ",\\\n";
                      } else {
                          ss << ")\n";
                      }
                  }
              } else {
                  // if we limit the number of colors, the new colors
                  // are the max_color appropriate interpolations in
                  // the original colormap
                  for (size_t i = 0; i < max_colors_; ++i) {
                      color_array c =
                          colormap_interpolation(static_cast<double>(i), 0., max_colors_ - 1.);
                      ss << "    " << i << "   " << c[0 + 1] << " " << c[1 + 1]
                         << " " << c[2 + 1];
                      if (i != max_colors_ - 1 || max_colors_ == 1) {
                          ss << ",\\\n";
                      } else {
                          ss << ")\n";
                      }
                  }
                  // we can't have just one color, at least in some versions
                  // of gnuplot we can't
                  if (max_colors_ == 1) {
                      color_array c =
                          colormap_interpolation(0, 0., max_colors_ - 1.);
                      ss << "    " << 1 << "   " << c[0 + 1] << " " << c[1 + 1]
                         << " " << c[2 + 1] << ")\n";
                  }
              }
              run_command(ss.str());
          }
      }
  
      std::tuple<double, double, double, double, double, double>
      axes_type::calculate_margins() {
          // There are conditions on which the xlim are smaller than requested
          // to make room for other elements
          double width_multiplier = is_3d() && !is_3d_map() ? 0.72 : 1.;
          double height_multiplier = is_3d() && !is_3d_map() ? 0.72 : 1.;
  
          double x = x_origin();
          double y = y_origin();
  
          // if parent has a title
          if (!parent_->title().empty()) {
              // open room for the title
              double room_needed = 0.2;
              double this_axes_contribution = room_needed * height();
              height_multiplier -= this_axes_contribution;
              // the further from y=0, the larger the shift
              y -= 1.5 * y * this_axes_contribution;
          }
  
          double rmargin_multiplier = (1. - width_multiplier) / 2.;
          double lmargin_multiplier = (1. - width_multiplier) / 2.;
          double tmargin_multiplier = (1. - height_multiplier) / 2.;
          double bmargin_multiplier = (1. - height_multiplier) / 2.;
  
          const auto bg_width = width_multiplier * width();
          const auto bg_height = height_multiplier * height();
  
          // We increase the margins for colorbar
          double colorbar_lmargin = 0.;
          double colorbar_rmargin = 0.;
          double colorbar_bmargin = 0.;
          double colorbar_tmargin = 0.;
  
          constexpr double extra_for_ticks = 0.05;
          if (cb_axis_.visible() && !cb_inside_) {
              // identify if colorbar is north/south/east/west
              auto [cbx, cby, cbw, cbh] = cb_position_;
              bool automatic_position = cbw == cbh;
              if (automatic_position) {
                  if (cb_vertical_) {
                      colorbar_rmargin = 0.04;
                  } else {
                      colorbar_bmargin = 0.1;
                  }
              } else if (cbw > cbh) {
                  // north/south
                  double south_margin = cby + cbh;
                  double north_margin = 1. - cby;
                  if (south_margin <= north_margin) {
                      // south
                      colorbar_bmargin = south_margin + extra_for_ticks;
                  } else {
                      // north
                      colorbar_tmargin = north_margin + extra_for_ticks;
                  }
              } else {
                  // west/east
                  double east_margin = cbx + cbw;
                  double west_margin = 1. - cbx;
                  if (east_margin <= west_margin) {
                      // east
                      colorbar_rmargin = east_margin + extra_for_ticks;
                  } else {
                      // west
                      colorbar_lmargin = west_margin + extra_for_ticks;
                  }
              }
          }
  
          const double lmargin =
              x + width() * lmargin_multiplier + colorbar_lmargin;
          const double rmargin =
              x + width() - width() * rmargin_multiplier - colorbar_rmargin;
          const double bmargin =
              y + height() * bmargin_multiplier + colorbar_bmargin;
          const double tmargin =
              y + height() - height() * tmargin_multiplier - colorbar_tmargin;
  
          return std::make_tuple(bg_width, bg_height,lmargin,rmargin,bmargin,tmargin);
      }
  
      void axes_type::run_position_margin_command() {
          include_comment("Axes position");
          run_command("set origin " + std::to_string(x_origin()) + "," +
                      std::to_string(y_origin()));
  
          auto [bg_width, bg_height, lmargin, rmargin, bmargin, tmargin] =
              calculate_margins();
  
          run_command("set size " + num2str(bg_width) + "," + num2str(bg_height));
  
          run_command("set lmargin at screen " + num2str(lmargin));
          run_command("set rmargin at screen " + num2str(rmargin));
          run_command("set bmargin at screen " + num2str(bmargin));
          run_command("set tmargin at screen " + num2str(tmargin));
  
          if (!axes_aspect_ratio_auto_) {
              if (is_3d()) {
                  run_command("set view equal xyz");
              } else {
                  run_command("set size ratio " +
                              std::to_string(axes_aspect_ratio_));
              }
          }
      }
  
      void axes_type::run_grid_command() {
          // if the user didn't explicitly set the grid, turn on the grid
          // by default in logscale
          if (!x_user_grid_) {
              x_grid_ = x_axis().scale() == axis_type::axis_scale::log || is_3d();
              x_minor_grid_ = x_axis().scale() == axis_type::axis_scale::log;
          }
          if (!y_user_grid_) {
              y_grid_ = y_axis().scale() == axis_type::axis_scale::log || is_3d();
              y_minor_grid_ = y_axis().scale() == axis_type::axis_scale::log;
          }
          if (!z_user_grid_) {
              z_grid_ = z_axis().scale() == axis_type::axis_scale::log || is_3d();
              z_minor_grid_ = z_axis().scale() == axis_type::axis_scale::log;
          }
          if (!r_user_grid_) {
              r_grid_ = r_axis().scale() == axis_type::axis_scale::log || is_polar();
              r_minor_grid_ = r_axis().scale() == axis_type::axis_scale::log;
          }
  
          // set grid xtics ytics mxtics mytics layerdefault linecolor 'blue'
          // linetype 1 linewidth 2, linecolor 'red'
          if (x_grid_ || x_minor_grid_ || y_grid_ || y_minor_grid_ || z_grid_ ||
              z_minor_grid_ || r_grid_ || r_minor_grid_) {
              include_comment("Create grid");
              std::string cmd = "set grid";
              if (x_grid_) {
                  cmd += " xtics";
              }
              if (x_minor_grid_) {
                  cmd += " mxtics";
              }
              if (y_grid_) {
                  cmd += " ytics";
              }
              if (y_minor_grid_) {
                  cmd += " mytics";
              }
              if (z_grid_) {
                  cmd += " ztics";
              }
              if (z_minor_grid_) {
                  cmd += " mztics";
              }
              if (r_grid_) {
                  cmd += " polar";
              }
              if (r_minor_grid_) {
                  cmd += " mrtics";
              }
              cmd += " vertical";
  
              if (grid_front_) {
                  cmd += " front";
              }
              cmd += grid_line_style_.plot_string(
                  line_spec::style_to_plot::plot_line_only, false);
              if (x_minor_grid_ || y_minor_grid_ || z_minor_grid_) {
                  cmd += ",";
                  if (grid_front_) {
                      cmd += " front ";
                  }
                  cmd += minor_grid_line_style_.plot_string(
                      line_spec::style_to_plot::plot_line_only, false);
              }
              run_command(cmd);
          }
      }
  
      void axes_type::run_box_command() {
          include_comment("Axes box");
          if (!visible_) {
              run_command("unset border");
              run_command("unset xtics");
              run_command("unset ytics");
              run_command("unset ztics");
              return;
          }
  
          /*  Border types for box
              1	bottom	    bottom left front
              2	left	    bottom left back
              4	top	        bottom right front
              8	right	    bottom right back
              16	no effect	left vertical
              32	no effect	back vertical
              64	no effect	right vertical
              128	no effect	front vertical
              256	no effect	top left back
              512	no effect	top right back
              1024	no effect	top left front
              2048	no effect	top right front
           */
          if (!box_ || is_polar()) {
              if (is_3d()) {
                  // create borders
                  // bottom left front 1
                  // bottom right front 4
                  // left vertical 16
                  run_command("set border 21 linecolor \"" +
                              to_string(x_axis().color()) + "\" linewidth " +
                              num2str(line_width()));
              } else if (is_2d()) {
                  // create left (2) and bottom (1)
                  unsigned border = 0;
                  if (x_axis().visible()) {
                      border += 1;
                  }
                  if (y_axis().visible()) {
                      border += 2;
                  }
                  run_command("set border " + num2str(border) + " linecolor \"" +
                              to_string(x_axis().color()) + "\" linewidth " +
                              num2str(line_width()));
              } else {
                  run_command("unset border");
              }
          } else {
              // if (box_)
              if (is_3d()) {
                  if (box_full_) {
                      run_command("set border 4095 linecolor \"" +
                                  to_string(box_color_) + "\" linewidth " +
                                  num2str(line_width()));
                  } else {
                      // create borders
                      // bottom left front 1
                      // bottom left back 2
                      // bottom right front 4
                      // bottom right back 8
                      // left vertical 16
                      // back vertical 32
                      // right vertical 64
                      // top left back 256
                      // top right back 512
                      run_command("set border 895 linecolor \"" +
                                  to_string(box_color_) + "\" linewidth " +
                                  num2str(line_width()));
                  }
              } else {
                  // create all borders (bottom / left / top / right) in black
                  run_command("set border 15 linecolor \"" +
                              to_string(box_color_) + "\" linewidth " +
                              num2str(line_width()));
              }
          }
      }
  
      void axes_type::run_axes_command() {
          // plot tics
          // http://www.gnuplot.info/docs_4.2/node182.html
          if (is_polar()) {
              run_command("set polar");
          }
          auto set_or_unset_axis = [this](class axis_type &ax, std::string axis_name,
                                          bool minor_ticks = false) {
              // cb is the only axis we don't unset if tics are empty
              if (ax.visible() && !ax.tick_values().empty()) {
                  if (ax.geographic()) {
                      run_command("set " + axis_name + "tics geographic");
                  }
                  if (ax.on_axis()) {
                      run_command("set " + axis_name + "tics axis");
                  }
                  run_command("set format " + axis_name + " \"" +
                              escape(ax.tick_label_format()) + "\"");
                  std::string tics_str = "set " + axis_name +
                                         "tics textcolor \"" +
                                         to_string(ax.color());
                  if (parent_->backend_->supports_fonts()) {
                      tics_str += "\" font \"" + font() + "," +
                                  num2str(unsigned(font_size())) + "\"";
                  }
                  tics_str +=
                      " nomirror in scale " + num2str(ax.tick_length()) + " ";
                  run_command(tics_str);
                  if (ax.tick_values_manual() || ax.tick_values().empty()) {
                      run_command("set " + axis_name + "tics" +
                                  ax.tick_values_string(minor_ticks) +
                                  ax.tick_rotate_string());
                  } else {
                      if (ax.tickangle() != 0.) {
                          run_command("set " + axis_name + "tics" +
                                      ax.tick_rotate_string());
                      }
                      if (minor_ticks) {
                          run_command("set m" + axis_name + "tics 2");
                      }
                  }
              } else {
                  if (axis_name == "r") {
                      // Hiding the r-axis works differently.
                      // We can't completely unset the tics or the grid does not
                      // appear properly
                      run_command("unset " + axis_name + "axis");
                      run_command("set " + axis_name + "tics scale 0");
                  } else if (axis_name == "cb") {
                      // the colorbar / colorbox has a special command to unset it
                      // if only ticks are empty we still want to show the box
                      if (!ax.visible()) {
                          run_command("unset colorbox");
                      }
                      if (ax.tick_values().empty()) {
                          run_command("unset cbtics");
                      }
  
                  } else {
                      run_command("unset " + axis_name + "tics");
                  }
              }
              const bool zero_axis_is_possible =
                  axis_name != "r" && axis_name != "t" && axis_name != "cb";
              if (zero_axis_is_possible) {
                  if (ax.zero_axis()) {
                      run_command("set " + axis_name + "zeroaxis linestyle -1");
                  } else {
                      run_command("unset " + axis_name + "zeroaxis");
                  }
              }
          };
  
          if (!visible_) {
              run_command("unset xtics");
              run_command("unset ytics");
              run_command("unset x2tics");
              run_command("unset y2tics");
              if (is_3d()) {
                  run_command("unset ztics");
              }
              if (is_polar()) {
                  run_command("unset raxis");
                  run_command("set rtics scale 0");
                  run_command("unset ttics");
              }
          } else {
              set_or_unset_axis(x_axis_, "x", x_minor_grid_);
              set_or_unset_axis(x2_axis_, "x2", x_minor_grid_);
              set_or_unset_axis(y_axis_, "y", y_minor_grid_);
              set_or_unset_axis(y2_axis_, "y2", y_minor_grid_);
  
              const bool is_map = azimuth_ == 0 && elevation_ == 90;
              if (is_3d() && !is_map) {
                  set_or_unset_axis(z_axis_, "z", z_minor_grid_);
              }
  
              if (is_polar()) {
                  set_or_unset_axis(r_axis_, "r", r_minor_grid_);
                  set_or_unset_axis(t_axis_, "t", t_minor_grid_);
              }
  
              if (cb_axis_.visible()) {
                  std::stringstream ss;
                  ss << "set colorbox";
                  if (!cb_axis_.reverse()) {
                      ss << " noinvert";
                  } else {
                      ss << " invert";
                  }
                  if (cb_vertical_) {
                      ss << " vertical";
                  } else {
                      ss << " horizontal";
                  }
                  if (cb_position_[2] == 0. && cb_position_[3] == 0.) {
                      ss << " default";
                  } else {
                      auto [x, y, w, h] = position_;
                      ss << " user";
                      ss << " origin " << x + cb_position_[0] * w << ","
                         << y + cb_position_[1] * h;
                      ss << " size " << cb_position_[2] * w << ","
                         << cb_position_[3] * h;
                  }
                  run_command(ss.str());
              }
              set_or_unset_axis(cb_axis_, "cb", false);
  
              /*
               * if (color_box_) {
                  run_command("set colorbox");
                  if (color_box_log_scale_) {
                      run_command("set logscale cb");
                  }
              } else {
                  run_command("unset colorbox");
                  if (color_box_log_scale_) {
                      run_command("unset logscale cb");
                  }
              }
               */
          }
  
          include_comment("Axes ranges");
          run_command("set xrange " + x_axis().range_string());
          run_command("set yrange " + y_axis().range_string());
          run_command("set y2range " + y2_axis().range_string());
          if (is_3d()) {
              run_command("set zrange " + z_axis().range_string());
          }
          if (is_polar() && r_axis().limits_mode_manual()) {
              run_command("set rrange " + r_axis().range_string());
          }
          if (cb_axis_.visible()) {
              run_command("set cbrange " + cb_axis_.range_string());
          }
  
          include_comment("Axes scale");
          if (x_axis().scale() == axis_type::axis_scale::log) {
              run_command("set logscale x");
          }
          if (y_axis().scale() == axis_type::axis_scale::log) {
              run_command("set logscale y");
          }
          if (is_3d()) {
              if (z_axis().scale() == axis_type::axis_scale::log) {
                  run_command("set logscale z");
              }
          }
  
          if (is_3d()) {
              include_comment("Axes view");
              run_command("set xyplane 0");
              // gnuplot azimuths are bounded [0:360]
              float gnuplot_azimuth_ = azimuth_ + 360;
              while (gnuplot_azimuth_ < 0) {
                  gnuplot_azimuth_ += 360;
              }
              while (gnuplot_azimuth_ > 360) {
                  gnuplot_azimuth_ -= 360;
              }
  
              // gnuplot elevation is bounded [0,180]
              // our elevation is bounded [-90,+90]
              // matlab -90 = gnuplot 180
              // matlab -45 = gnuplot 135
              // matlab  0  = gnuplot  90
              // matlab +45 = gnuplot  45
              // matlab +90 = gnuplot   0
              float gnuplot_elevation_ = 90 - elevation_;
              while (gnuplot_elevation_ < 0) {
                  gnuplot_elevation_ += 180;
              }
              while (gnuplot_elevation_ > 180) {
                  gnuplot_elevation_ -= 180;
              }
              run_command("set view " + num2str(gnuplot_elevation_) + "," +
                          num2str(gnuplot_azimuth_));
          }
      }
  
      void axes_type::run_labels_command() {
          include_comment("Axes labels");
          //      set xlabel {"<label>"} {offset <offset>} {font
          //      "<font>{,<size>}"}
          //                {{textcolor | tc} {lt <line_type> | default}}
          //                {{no}enhanced} {rotate by <degrees>}
          if (!x_axis_.label().empty()) {
              std::string cmd = "set xlabel" + x_axis_.label_string();
              run_command(cmd);
          }
          if (!y_axis_.label().empty()) {
              std::string cmd = "set ylabel" + y_axis_.label_string();
              run_command(cmd);
          }
          if (!y2_axis_.label().empty()) {
              std::string cmd = "set y2label" + y2_axis_.label_string();
              run_command(cmd);
          }
  
          const bool is_map = azimuth_ == 0 && elevation_ == 90;
          if (is_3d() && !z_axis_.label().empty() && !is_map) {
              std::string cmd =
                  "set zlabel" + z_axis_.label_string() + " rotate parallel";
              run_command(cmd);
          }
          if (!cb_axis_.label().empty()) {
              std::string cmd =
                  "set cblabel" + cb_axis_.label_string() + " rotate parallel";
              run_command(cmd);
          }
      }
  
      void axes_type::run_legend_command() {
          //  set key {on|off} {default}
          //  {{inside | outside} | {lmargin | rmargin | tmargin | bmargin}
          //      | {at <position>}}
          //  {left | right | center} {top | bottom | center}
          //  {vertical | horizontal} {Left | Right}
          //  {{no}reverse} {{no}invert}
          //  {samplen <sample_length>} {spacing <vertical_spacing>}
          //  {width <width_increment>}
          //  {height <height_increment>}
          //  {{no}autotitle {columnheader}}
          //  {title "<text>"} {{no}enhanced}
          //  {{no}box { {linestyle | ls <line_style>}
          //              | {linetype | lt <line_type>}
          //              {linewidth | lw <line_width>}}}
          include_comment("Axes legend");
          // Gnuplot version needs to be 5.2.6+ for keyentry
          bool ok = true;
          if (parent_->backend_->consumes_gnuplot_commands()) {
              std::pair<int, int> v = backend::gnuplot::gnuplot_version();
              if (v.first < 5 || v.second < 2) {
                  ok = false;
              }
          }
          if (legend_ == nullptr || !legend_->visible() || !ok) {
              run_command("set key off");
          } else {
              bool legends_to_plot = !legend_->empty();
              if (!legends_to_plot) {
                  for (const auto &child : children_) {
                      if (!child->display_name().empty()) {
                          legends_to_plot = true;
                          break;
                      }
                  }
              }
  
              if (!legends_to_plot) {
                  run_command("set key off");
              } else {
                  std::string cmd = "set key on";
                  if (parent_->backend_->supports_fonts()) {
                      cmd += " font \"" + legend_->font_name() + "," +
                             num2str(legend_->font_size()) +
                             "\" textcolor rgb '" +
                             to_string(legend_->text_color()) + "'";
                  }
                  // position
                  cmd += legend_->inside() ? " inside" : " outside";
                  if (legend_->manual_position()) {
                      cmd += " at " + num2str(legend_->position()[0]) + "," +
                             num2str(legend_->position()[1]);
                  }
                  switch (legend_->horizontal_location()) {
                  case legend::horizontal_alignment::left:
                      cmd += " left";
                      break;
                  case legend::horizontal_alignment::right:
                      cmd += " right";
                      break;
                  case legend::horizontal_alignment::center:
                      cmd += " center";
                      break;
                  }
                  switch (legend_->vertical_location()) {
                  case legend::vertical_alignment::top:
                      cmd += " top";
                      break;
                  case legend::vertical_alignment::bottom:
                      cmd += " bottom";
                      break;
                  case legend::vertical_alignment::center:
                      cmd += " center";
                      break;
                  }
                  cmd += legend_->vertical() ? " vertical" : " horizontal";
                  // text aligned to the left
                  cmd += " Left";
                  cmd += " opaque";
                  cmd +=
                      legend_->label_after_sample() ? " reverse" : " noreverse";
                  cmd += legend_->invert() ? " invert" : " noinvert";
                  // increase the width to make labels fit
                  cmd += " width 1";
                  // forbid automatic labels
                  // we are going to generate them explicitly for our xlim objects
                  cmd += " noautotitle";
                  // box
                  if (!legend_->title().empty()) {
                      cmd += " title \"";
                      if (iequals(legend_->font_weight(), "bold")) {
                          cmd += "{/:Bold ";
                      }
                      cmd += escape(legend_->title());
                      if (iequals(legend_->font_weight(), "bold")) {
                          cmd += "}";
                      }
                      cmd += "\"";
                      if (parent_->backend_->supports_fonts()) {
                          cmd += " font \"" + escape(legend_->font_name()) + "," +
                                 num2str(unsigned(legend_->font_size())) + "\"";
                          cmd += " textcolor rgb \"" +
                                 to_string(legend_->text_color()) + "\"";
                      }
                  }
                  if (legend_->box() && legend_->box_line().has_line()) {
                      cmd += " box";
                      cmd += legend_->box_line().plot_string(
                          line_spec::style_to_plot::plot_line_only, false);
                  } else {
                      cmd += " nobox";
                  }
                  if (legend_->num_columns() != 0) {
                      size_t n_strings = legend_->size();
                      for (size_t i = 0; i < children_.size(); ++i) {
                          n_strings += !children_[i]->display_name().empty();
                      }
                      cmd += " maxrows " +
                             num2str(ceil(static_cast<double>(n_strings) /
                                          legend_->num_columns()));
                  }
                  if (legend_->num_rows() != 0) {
                      cmd += " maxrows " + num2str(legend_->num_rows());
                  }
  
                  run_command(cmd);
              }
          }
      }
  
      void axes_type::run_background_command() {
          const bool saving_to_file = !parent_->backend_->output().empty();
          const bool using_default_color =
              color_ == std::array<float, 4>{0, 1.0, 1.0, 1.0};
          if (saving_to_file && using_default_color) {
              return;
          }
          include_comment("Axes background");
          // plot background
          if (visible_) {
              // rectangle behind the graph
              if (color_ != parent()->color()) {
                  if (is_2d() || is_3d_map()) {
                      run_command("set object 2 rectangle from graph 0,0 to "
                                  "graph 1,1 behind fillcolor rgb \"" +
                                  to_string(color()) +
                                  "\" fillstyle solid 1.0 noborder");
                  } else if (is_3d()) {
                      auto v = backend::gnuplot::gnuplot_version();
                      const bool has_wall_option =
                          v.first > 5 || (v.first == 5 && v.second >= 5);
                      if (has_wall_option) {
                          run_command("set wall z0 fs transparent solid 0.5 "
                                      "border -1 fc 'slategray'");
                          run_command("set wall x0 fs transparent solid 0.5 "
                                      "border -1 fc 'slategray'");
                          run_command("set wall y1 fs transparent solid 0.5 "
                                      "border -1 fc 'slategray'");
                      }
                  } else if (is_polar()) {
                      double m = children_[0]->xmax();
                      run_command("set object 2 ellipse at first 0,0 size " +
                                  num2str(m * 2) + "," + num2str(m * 2) +
                                  " angle 0 behind fillcolor rgb \"" +
                                  to_string(color()) +
                                  "\" fillstyle solid 1.0 noborder");
                  }
              }
              // rectangle behind the whole area allocated for the xlim
              // we don't use it anymore
              /*
               * if (false && color_ != parent()->color()) {
               *  run_command("set object 2 rectangle from screen " +
               * std::to_string(x_origin()) + "," + std::to_string(y_origin()) + "
               * to screen " + std::to_string(x_origin() + width()) + "," +
               *              std::to_string(y_origin() + height()) + " behind
               * fillcolor rgb \"" + to_string(color_outside_) + "\" fillstyle
               * solid 1.0 noborder");
               * }
               */
          }
      }
  
      void axes_type::run_empty_plot_command() {
          include_comment("Empty xlim");
          if (x_axis().limits_mode_auto()) {
              run_command("set yrange [0:1]");
          }
          if (y_axis().limits_mode_auto()) {
              run_command("set yrange [0:1]");
          }
          run_command("set key off");
          if (y_axis().limits_mode_auto() ||
              !std::isfinite(y_axis().limits()[1])) {
              run_command("plot 2 with lines");
          } else {
              run_command("plot " + std::to_string(y_axis().limits_[1] + 1) +
                          " with lines");
          }
      }
  
      void axes_type::run_plot_objects_command() {
          include_comment("Axes objects");
          // Plot empty xlim if there are no objects
          if (children_.empty()) {
              run_empty_plot_command();
              return;
          }
  
          // Set variable commands
          std::string set_variables_command;
          for (const auto &child : children_) {
              set_variables_command += child->set_variables_string();
          }
          run_command(set_variables_command);
  
          // Plot command
          std::string plot_command = !is_3d() ? "plot " : "splot ";
  
          // Plot all children
          bool first = true;
          for (const auto &child : children_) {
              if (!first) {
                  plot_command += ",\\\n         ";
              } else {
                  first = false;
              }
              plot_command += child->plot_string();
          }
  
          // Keyentry commands (legends for each child)
          if (legend_ != nullptr && legend_->visible()) {
              bool ok = true;
              static bool msg_shown_once = false;
              // Gnuplot version needs to be 5.2.6+ for keyentry
              if (parent_->backend_->consumes_gnuplot_commands()) {
                  std::pair<int, int> v = backend::gnuplot::gnuplot_version();
                  if (v.first < 5 || v.second < 2) {
                      if (!msg_shown_once) {
                          std::cerr
                              << "You need gnuplot 5.2.6+ to include legends"
                              << std::endl;
                          msg_shown_once = true;
                      }
                      ok = false;
                  }
              }
  
              if (ok) {
                  auto legend_it = legend_->begin();
                  auto legend_end = legend_->end();
                  for (auto child_it = children_.begin();
                       child_it != children_.end(); ++child_it) {
                      const bool no_legend = legend_it == legend_end &&
                                             (*child_it)->display_name().empty();
                      if (no_legend) {
                          continue;
                      }
                      const auto &child = *child_it;
                      if (!first) {
                          plot_command += ", ";
                      } else {
                          first = false;
                      }
                      // https://stackoverflow.com/a/60624922/2983585
                      // http://gnuplot.sourceforge.net/demo_svg_5.5/custom_key.html
                      plot_command += child->legend_string(legend_it, legend_end);
                  }
              }
          }
  
          if (plot_command != "plot ") {
              run_command(plot_command);
          } else {
              run_empty_plot_command();
          }
  
          // data commands
          for (const auto &child : children_) {
              run_command(child->data_string());
          }
  
          // unset variables command
          std::string unset_variables_command;
          for (const auto &child : children_) {
              unset_variables_command += child->unset_variables_string();
          }
          run_command(unset_variables_command);
      }
  
      void axes_type::run_commands() {
          run_colormap_command();
          run_position_margin_command();
          run_title_command();
          run_box_command();
          run_grid_command();
          run_axes_command();
          run_labels_command();
          run_legend_command();
          run_background_command();
          run_plot_objects_command();
      }
  
      void axes_type::run_draw_commands() {
          run_background_draw_commands();
          run_grid_draw_commands();
          run_box_draw_commands();
          run_axes_draw_commands();
          run_labels_draw_commands();
          run_title_draw_commands();
          run_legend_draw_commands();
          run_plot_objects_draw_commands();
      }
  
      void axes_type::run_background_draw_commands() {
          auto [w, h, lm, rm, bm, tm] = calculate_margins();
          double view_width = parent_->backend_->width();
          double x1 = lm * view_width;
          double x2 = rm * view_width;
          double view_height = parent_->backend_->height();
          double y1 = bm * view_height;
          double y2 = tm * view_height;
          parent_->backend_->draw_rectangle(x1, x2, y1, y2, this->color_);
      }
  
      void axes_type::run_title_draw_commands() {
          // This should draw the title on top of the axes
      }
  
      void axes_type::run_box_draw_commands() {
          auto [w, h, lm, rm, bm, tm] = calculate_margins();
          double view_width = parent_->backend_->width();
          double x1 = lm * view_width;
          double x2 = rm * view_width;
          double view_height = parent_->backend_->height();
          double y1 = bm * view_height;
          double y2 = tm * view_height;
          const std::array<float, 4> color = {0., 0., 0., 0.};
          std::vector<double> box_xs = {x1, x2, x2, x1, x1};
          std::vector<double> box_ys = {y1, y1, y2, y2, y1};
          parent_->backend_->draw_path(box_xs, box_ys, color);
      }
  
      void axes_type::run_grid_draw_commands() {}
  
      void axes_type::run_axes_draw_commands() {
          auto [w, h, lm, rm, bm, tm] = calculate_margins();
          double view_width = parent_->backend_->width();
          double view_height = parent_->backend_->height();
          // to draw the x axis we
          // - normalize x values
          // - keep y as absolute values
          auto xlimits = xlim();
          auto ylimits = ylim();
          std::vector<double> cx;
          if (x_axis_.tick_values_manual()) {
              cx = x_axis_.tick_values();
          } else {
              ticks_results xticks_results = calcticks(xlimits[0], xlimits[1], true, 1.25, true, false);
              cx = xticks_results.ticks;
          }
          // clamp values outside x limits
          cx = transform(cx, [&](double x) { return std::clamp(x,xlimits[0],xlimits[1]); });
  
          // convert x values to viewport range
          double viewport_xmin = lm * view_width;
          double viewport_xmax = rm * view_width;
          double viewport_xrange = viewport_xmax - viewport_xmin;
  
          auto xrange = xlimits[1] - xlimits[0];
          for (auto &v : cx) {
              v -= xlimits[0];
              v /= xrange;
              v *= viewport_xrange;
              v += viewport_xmin;
          }
  
          // get viewport range for y axis
          double viewport_ymin = bm * view_height;
          double viewport_ymax = tm * view_height;
          double viewport_yrange = viewport_ymax - viewport_ymin;
  
          // The meaning of tick length is proportional
          // to the largest of x and y viewport size.
          // This unit makes the tick sizes look the same on x and y
          // while still proportional to the axes object.
          double tick_length_multiplier = std::max(viewport_xrange, viewport_yrange) * 0.015;
  
          // draw a path for each x tick
          double xtick_length_multiplier = std::min(tick_length_multiplier, viewport_xrange);
          for (auto &v : cx) {
              parent_->backend_->draw_path({v,v}, {viewport_ymin, viewport_ymin + x_axis_.tick_length() * xtick_length_multiplier},x_axis_.color_);
              parent_->backend_->draw_path({v,v}, {viewport_ymax, viewport_ymax - x_axis_.tick_length() * xtick_length_multiplier},x_axis_.color_);
          }
  
          // to draw the y axis we
          // - normalize y values
          // - keep x as absolute values
          std::vector<double> cy;
          if (y_axis_.tick_values_manual()) {
              cy = y_axis_.tick_values();
          } else {
              ticks_results yticks_results = calcticks(ylimits[0], ylimits[1], false, 1.25, true, false);
              cy = yticks_results.ticks;
          }
          // clamp values outside x limits
          cy = transform(cy, [&](double x) { return std::clamp(x,ylimits[0],ylimits[1]); });
  
          // convert y values to viewport range
          auto yrange = ylimits[1] - ylimits[0];
          for (auto &v : cy) {
              v -= ylimits[0];
              v /= yrange;
              v *= viewport_yrange;
              v += viewport_ymin;
          }
  
          // draw a path for each y tick
          double ytick_length_multiplier = std::min(tick_length_multiplier, viewport_yrange);
          for (auto &v : cy) {
              parent_->backend_->draw_path({viewport_xmin, viewport_xmin + y_axis_.tick_length() * ytick_length_multiplier}, {v,v}, y_axis_.color_);
              parent_->backend_->draw_path({viewport_xmax, viewport_xmax - y_axis_.tick_length() * ytick_length_multiplier}, {v,v}, y_axis_.color_);
          }
      }
  
      void axes_type::run_labels_draw_commands() {}
  
      void axes_type::run_legend_draw_commands() {}
  
      void axes_type::run_plot_objects_draw_commands() {
          for (const auto &child : children_) {
              child->run_draw_commands();
          }
      }
  
      void axes_type::run_command(const std::string &command) {
          parent_->run_command("    " + command);
      }
  
      void axes_type::include_comment(const std::string &command) {
          parent_->run_command("    # " + command);
      }
  
      bool axes_type::is_3d() {
          return children_.empty() ? false : children_[0]->is_3d();
      }
  
      bool axes_type::is_3d_map() {
          return children_.empty() ? false : children_[0]->is_3d_map();
      }
  
      bool axes_type::is_2d() {
          return children_.empty() ? true : children_[0]->is_2d();
      }
  
      bool axes_type::is_polar() {
          return children_.empty() ? false : children_[0]->is_polar();
      }
  
      axes_handle axes_type::copy(figure_handle parent) {
          axes_handle new_axes = std::make_shared<class axes_type>(*this);
          new_axes->parent_ = parent.get();
          new_axes->parent_->add_axes(new_axes, false, false);
          return new_axes;
      }
  
      std::array<double, 4> axes_type::child_limits() const {
          if (children_.empty()) {
              return {-10, +10, -10, +10};
          }
          auto xminit =
              std::min_element(children_.begin(), children_.end(),
                               [](axes_object_handle a, axes_object_handle b) {
                                   return a->xmin() < b->xmin();
                               });
          auto xmaxit =
              std::max_element(children_.begin(), children_.end(),
                               [](axes_object_handle a, axes_object_handle b) {
                                   return a->xmax() < b->xmax();
                               });
          auto yminit =
              std::min_element(children_.begin(), children_.end(),
                               [](axes_object_handle a, axes_object_handle b) {
                                   return a->ymin() < b->ymin();
                               });
          auto ymaxit =
              std::max_element(children_.begin(), children_.end(),
                               [](axes_object_handle a, axes_object_handle b) {
                                   return a->ymax() < b->ymax();
                               });
          return {(*xminit)->xmin(), (*xmaxit)->xmax(), (*yminit)->ymin(),
                  (*ymaxit)->ymax()};
      }
  
      axes_type::axes_type() : axes_type(nullptr) {
          t_axis_.tick_label_format_ = "%gº";
          t_axis_.tick_values_ = iota(0, 30, 330);
          t_axis_.tick_values_automatic_ = false;
      }
  
      axes_type::axes_type(class figure_type *parent)
          : axes_type(parent, default_axes_position) {
          t_axis_.tick_label_format_ = "%gº";
          t_axis_.tick_values_ = iota(0, 30, 330);
          t_axis_.tick_values_automatic_ = false;
          font_ = parent->font();
      }
  
      axes_type::axes_type(class figure_type *parent,
                           std::array<float, 4> position)
          : parent_(parent), position_(position), x_axis_(this, -10, +10, true),
            x2_axis_(this, inf, inf, false), y_axis_(this, inf, inf, true),
            y2_axis_(this, inf, inf, false), z_axis_(this, inf, inf, true) {
          grid_line_style_.color({0.85, 0.15, 0.15, 0.15});
          minor_grid_line_style_.color({0.9, 0.1, 0.1, 0.1});
          t_axis_.tick_label_format_ = "%gº";
          t_axis_.tick_values_ = iota(0, 30, 330);
          t_axis_.tick_values_automatic_ = false;
          font_ = parent->font();
      }
  
      axes_type::axes_type(figure_handle parent) : axes_type(parent.get()) {
          t_axis_.tick_label_format_ = "%gº";
          t_axis_.tick_values_ = iota(0, 30, 330);
          t_axis_.tick_values_automatic_ = false;
          font_ = parent->font();
      }
  
      axes_type::axes_type(figure_handle parent, std::array<float, 4> position)
          : axes_type(parent.get(), position) {
          t_axis_.tick_label_format_ = "%gº";
          t_axis_.tick_values_ = iota(0, 30, 330);
          t_axis_.tick_values_automatic_ = false;
          font_ = parent->font();
      }
  
      const class axis_type &axes_type::x_axis() const { return x_axis_; }
  
      class axis_type &axes_type::x_axis() {
          return x_axis_;
      }
  
      const class axis_type &axes_type::x2_axis() const { return x2_axis_; }
  
      class axis_type &axes_type::x2_axis() {
          return x2_axis_;
      }
  
      const class axis_type &axes_type::y_axis() const { return y_axis_; }
  
      const std::string &axes_type::xlabel() const { return x_axis_.label(); }
  
      void axes_type::xlabel(std::string_view str) {
          x_axis_.label(str);
          touch();
      }
  
      const std::string &axes_type::x2label() const { return x2_axis_.label(); }
  
      void axes_type::x2label(std::string_view str) {
          x2_axis_.label(str);
          touch();
      }
  
      const std::string &axes_type::ylabel() const { return y_axis_.label(); }
  
      void axes_type::ylabel(std::string_view str) {
          y_axis_.label(str);
          touch();
      }
  
      const std::string &axes_type::xtickformat() const {
          return x_axis_.tick_label_format();
      }
  
      void axes_type::xtickformat(std::string_view str) {
          x_axis_.tick_label_format(str);
      }
  
      const std::string &axes_type::x2tickformat() const {
          return x2_axis_.tick_label_format();
      }
  
      void axes_type::x2tickformat(std::string_view str) {
          x2_axis_.tick_label_format(str);
      }
  
      const std::string &axes_type::ytickformat() const {
          return y_axis_.tick_label_format();
      }
  
      void axes_type::ytickformat(std::string_view str) {
          y_axis_.tick_label_format(str);
      }
  
      const std::string &axes_type::y2tickformat() const {
          return y2_axis_.tick_label_format();
      }
  
      void axes_type::y2tickformat(std::string_view str) {
          y2_axis_.tick_label_format(str);
      }
  
      const std::string &axes_type::ztickformat() const {
          return z_axis_.tick_label_format();
      }
  
      void axes_type::ztickformat(std::string_view str) {
          z_axis_.tick_label_format(str);
      }
  
      const std::string &axes_type::rtickformat() const {
          return r_axis_.tick_label_format();
      }
  
      void axes_type::rtickformat(std::string_view str) {
          r_axis_.tick_label_format(str);
      }
  
      const std::string &axes_type::y2label() const { return y2_axis_.label(); }
  
      void axes_type::y2label(std::string_view str) {
          y2_axis_.label(str);
          touch();
      }
  
      const std::string &axes_type::zlabel() const { return z_axis_.label(); }
  
      void axes_type::zlabel(std::string_view str) {
          z_axis_.label(str);
          touch();
      }
  
      const std::string &axes_type::rlabel() const { return r_axis_.label(); }
  
      void axes_type::rlabel(std::string_view str) {
          r_axis_.label(str);
          touch();
      }
  
      const std::vector<double> &axes_type::xticks() const {
          return x_axis_.tick_values();
      }
  
      void axes_type::xticks(const std::vector<double> &ticks) {
          x_axis_.tick_values(ticks);
      }
  
      const std::vector<double> &axes_type::x2ticks() const {
          return x2_axis_.tick_values();
      }
  
      void axes_type::x2ticks(const std::vector<double> &ticks) {
          x2_axis_.tick_values(ticks);
      }
  
      const std::vector<double> &axes_type::yticks() const {
          return y_axis_.tick_values();
      }
  
      void axes_type::yticks(const std::vector<double> &ticks) {
          y_axis_.tick_values(ticks);
      }
  
      const std::vector<double> &axes_type::y2ticks() const {
          return y2_axis_.tick_values();
      }
  
      void axes_type::y2ticks(const std::vector<double> &ticks) {
          y2_axis_.tick_values(ticks);
      }
  
      const std::vector<double> &axes_type::zticks() const {
          return z_axis_.tick_values();
      }
  
      void axes_type::zticks(const std::vector<double> &ticks) {
          z_axis_.tick_values(ticks);
      }
  
      const std::vector<double> &axes_type::rticks() const {
          return r_axis_.tick_values();
      }
  
      void axes_type::rticks(const std::vector<double> &ticks) {
          r_axis_.tick_values(ticks);
      }
  
      class axis_type &axes_type::y_axis() {
          return y_axis_;
      }
  
      const class axis_type &axes_type::y2_axis() const { return y2_axis_; }
  
      class axis_type &axes_type::y2_axis() {
          return y2_axis_;
      }
  
      const class axis_type &axes_type::z_axis() const { return z_axis_; }
  
      class axis_type &axes_type::z_axis() {
          return z_axis_;
      }
  
      const class axis_type &axes_type::r_axis() const { return r_axis_; }
  
      class axis_type &axes_type::r_axis() {
          return r_axis_;
      }
  
      const std::string &axes_type::font() const {
          if (font_) {
              return *font_;
          } else {
              return parent_->font();
          }
      }
  
      void axes_type::font(std::string_view font) {
          font_ = font;
          touch();
      }
  
      float axes_type::font_size() const {
          if (font_size_) {
              return *font_size_;
          } else {
              return parent_->font_size();
          }
      }
  
      void axes_type::font_size(float font_size) {
          font_size_ = font_size;
          touch();
      }
  
      const std::string &axes_type::font_weight() const { return font_weight_; }
  
      void axes_type::font_weight(std::string_view font_weight) {
          font_weight_ = font_weight;
          touch();
      }
  
      float axes_type::title_font_size_multiplier() const {
          return title_font_size_multiplier_;
      }
  
      void
      axes_type::title_font_size_multiplier(float title_font_size_multiplier) {
          title_font_size_multiplier_ = title_font_size_multiplier;
          touch();
      }
  
      const std::string &axes_type::title_font_weight() const {
          return title_font_weight_;
      }
  
      void axes_type::title_font_weight(std::string_view title_font_weight) {
          title_font_weight_ = title_font_weight;
          touch();
      }
  
      const color_array &axes_type::color() const { return color_; }
  
      void axes_type::color(const color_array &color) {
          color_ = color;
          touch();
      }
  
      void axes_type::color(const std::array<float, 3> &color) {
          color_ = {0, color[0], color[1], color[2]};
          touch();
      }
  
      void axes_type::color(std::string_view c) { color(string_to_color(c)); }
  
      void axes_type::color(const enum color &c) { color(to_array(c)); }
  
      const std::vector<color_array> &axes_type::colororder() const {
          return colororder_;
      }
  
      void axes_type::colororder(const std::vector<color_array> &colororder) {
          colororder_ = colororder;
          touch();
      }
  
      size_t axes_type::colororder_index() const { return colororder_index_; }
  
      void axes_type::colororder_index(size_t colororder_index) {
          colororder_index_ = colororder_index;
      }
  
      void axes_type::bump_colororder_index() {
          colororder_index_ = (colororder_index_ + 1) % colororder_.size();
      }
  
      color_array axes_type::get_color_and_bump() {
          color_array r = colororder_[colororder_index_];
          bump_colororder_index();
          return r;
      }
  
      const color_array axes_type::grid_color() const {
          return grid_line_style_.color();
      }
  
      void axes_type::grid_color(const color_array &grid_color) {
          grid_line_style_.color(grid_color);
          touch();
      }
  
      float axes_type::grid_alpha() const { return grid_line_style_.alpha(); }
  
      void axes_type::grid_alpha(float grid_alpha) {
          grid_line_style_.alpha(grid_alpha);
          touch();
      }
  
      bool axes_type::handle_visibility() const { return handle_visibility_; }
  
      void axes_type::handle_visibility(bool handle_visibility) {
          handle_visibility_ = handle_visibility;
          touch();
      }
  
      float axes_type::minor_grid_alpha() const {
          return minor_grid_line_style_.alpha();
      }
  
      void axes_type::minor_grid_alpha(float minor_grid_alpha) {
          minor_grid_line_style_.alpha(minor_grid_alpha);
          touch();
      }
  
      const color_array axes_type::minor_grid_color() const {
          return minor_grid_line_style_.color();
      }
  
      void axes_type::minor_grid_color(const color_array &minor_grid_color) {
          minor_grid_line_style_.color(minor_grid_color);
          touch();
      }
  
      line_spec axes_type::minor_grid_line_style() const {
          return minor_grid_line_style_;
      }
  
      void axes_type::minor_grid_line_style(line_spec minor_grid_line_style) {
          minor_grid_line_style_ = minor_grid_line_style;
          touch();
      }
  
      line_spec axes_type::grid_line_style() const { return grid_line_style_; }
  
      void axes_type::grid_line_style(line_spec grid_line_style) {
          grid_line_style_ = grid_line_style;
          touch();
      }
  
      const std::string &axes_type::line_style_order() const {
          return line_style_order_;
      }
  
      void axes_type::line_style_order(std::string_view line_style_order) {
          line_style_order_ = line_style_order;
          touch();
      }
  
      size_t axes_type::line_style_order_index() const {
          return line_style_order_index_;
      }
  
      void axes_type::line_style_order_index(size_t line_style_order_index) {
          line_style_order_index_ = line_style_order_index;
          touch();
      }
  
      float axes_type::line_width() const { return line_width_; }
  
      void axes_type::line_width(float line_width) {
          line_width_ = line_width;
          touch();
      }
  
      bool axes_type::next_plot_replace() const { return next_plot_replace_; }
  
      void axes_type::next_plot_replace(bool next_plot_replace) {
          if (next_plot_replace != next_plot_replace_) {
              next_plot_replace_ = next_plot_replace;
          }
      }
  
      bool axes_type::hold() const { return !next_plot_replace(); }
  
      void axes_type::hold(bool v) { next_plot_replace(!v); }
  
      const std::vector<std::shared_ptr<class axes_object>> &
      axes_type::children() const {
          return children_;
      }
  
      void axes_type::children(
          const std::vector<std::shared_ptr<axes_object>> &children) {
          children_ = children;
          touch();
      }
  
      const legend_handle &axes_type::legend() const { return legend_; }
  
      void axes_type::legend(const legend_handle &legend) {
          legend_ = legend;
          touch();
      }
  
      void axes_type::legend(const std::vector<std::string> &names) {
          if (legend_ == nullptr) {
              legend_.reset(new class legend(this, names));
          } else {
              legend_->strings(names);
          }
          touch();
      }
  
      const class figure_type *axes_type::parent() const { return parent_; }
  
      class figure_type *axes_type::parent() {
          return parent_;
      }
  
      void axes_type::parent(class figure_type *p) {
          parent_ = p;
          touch();
      }
  
      const std::string &axes_type::title() const { return title_; }
  
      void axes_type::title(std::string_view title) {
          title_ = title;
          title_visible_ = !title.empty();
          if (parent_->children_.size() == 1 && !parent_->quiet_mode()) {
              run_title_command();
              run_command("replot");
              parent_->flush_commands();
          } else {
              touch();
          }
      }
  
      const std::array<float, 4> &axes_type::position() const {
          return position_;
      }
  
      void axes_type::position(const std::array<float, 4> &position) {
          position_ = position;
          touch();
      }
  
      const float &axes_type::x_origin() const { return position_[0]; }
  
      void axes_type::x_origin(float x) {
          position_[0] = x;
          touch();
      }
  
      const float &axes_type::y_origin() const { return position_[1]; }
  
      void axes_type::y_origin(float y) {
          position_[1] = y;
          touch();
      }
  
      const float &axes_type::width() const { return position_[2]; }
  
      void axes_type::width(float w) {
          position_[2] = w;
          touch();
      }
  
      const float &axes_type::height() const { return position_[3]; }
  
      void axes_type::height(float h) {
          position_[3] = h;
          touch();
      }
  
      bool axes_type::box() const { return box_; }
  
      void axes_type::box(bool box) {
          if (box != box_) {
              box_ = box;
              touch();
          }
      }
  
      bool axes_type::box_full() const { return box_full_; }
  
      void axes_type::box_full(bool box_full) {
          if (box_full != box_full_) {
              box_full_ = box_full;
              if (box_full) {
                  box_ = true;
              }
              touch();
          }
      }
  
      bool axes_type::title_visible() const { return title_visible_; }
  
      void axes_type::title_visible(bool title_visible) {
          title_visible_ = title_visible;
          touch();
      }
  
      bool axes_type::visible() const { return visible_; }
  
      void axes_type::visible(bool visible) {
          visible_ = visible;
          touch();
      }
  
      void axes_type::limits_mode(keyword_automatic_type automatic) {
          limits_mode_automatic(true);
          touch();
      }
  
      void axes_type::limits_mode(keyword_manual_type manual) {
          limits_mode_automatic(false);
          touch();
      }
  
      void axes_type::limits_mode_automatic(bool automatic) {
          x_axis_.limits_mode_auto(automatic);
          y_axis_.limits_mode_auto(automatic);
          y2_axis_.limits_mode_auto(automatic);
          z_axis_.limits_mode_auto(automatic);
          touch();
      }
  
      bool axes_type::limits_mode_automatic() const {
          return x_axis_.limits_mode_auto() && y_axis_.limits_mode_auto() &&
                 y2_axis_.limits_mode_auto() && z_axis_.limits_mode_auto();
      }
  
      void axes_type::limits_mode_manual(bool manual) {
          x_axis_.limits_mode_auto(!manual);
          y_axis_.limits_mode_auto(!manual);
          y2_axis_.limits_mode_auto(!manual);
          z_axis_.limits_mode_auto(!manual);
          touch();
      }
  
      bool axes_type::limits_mode_manual() const {
          return x_axis_.limits_mode_manual() && y_axis_.limits_mode_manual() &&
                 y2_axis_.limits_mode_manual() && z_axis_.limits_mode_manual();
      }
  
      bool axes_type::grid() const { return x_grid_ && y_grid_ && z_grid_; }
  
      void axes_type::grid(bool grid) {
          x_grid_ = grid;
          x_user_grid_ = true;
          y_grid_ = grid;
          y_user_grid_ = true;
          z_grid_ = grid;
          z_user_grid_ = true;
          touch();
      }
  
      bool axes_type::minor_grid() const {
          return x_minor_grid_ && y_minor_grid_ && z_minor_grid_;
      }
  
      void axes_type::minor_grid(bool minor_grid) {
          x_minor_grid_ = minor_grid;
          x_user_grid_ = true;
          y_minor_grid_ = minor_grid;
          y_user_grid_ = true;
          z_minor_grid_ = minor_grid;
          z_user_grid_ = true;
          touch();
      }
  
      bool axes_type::x_grid() const { return x_grid_; }
  
      void axes_type::x_grid(bool x_grid) {
          x_grid_ = x_grid;
          x_user_grid_ = true;
          touch();
      }
  
      bool axes_type::x_minor_grid() const { return x_minor_grid_; }
  
      void axes_type::x_minor_grid(bool x_minor_grid) {
          x_minor_grid_ = x_minor_grid;
      }
  
      bool axes_type::y_grid() const { return y_grid_; }
  
      void axes_type::y_grid(bool y_grid) {
          y_grid_ = y_grid;
          x_user_grid_ = true;
          touch();
      }
  
      bool axes_type::y_minor_grid() const { return y_minor_grid_; }
  
      void axes_type::y_minor_grid(bool y_minor_grid) {
          y_minor_grid_ = y_minor_grid;
      }
  
      bool axes_type::z_grid() const { return z_grid_; }
  
      void axes_type::z_grid(bool z_grid) {
          z_grid_ = z_grid;
          x_user_grid_ = true;
          touch();
      }
  
      bool axes_type::z_minor_grid() const { return z_minor_grid_; }
  
      void axes_type::z_minor_grid(bool z_minor_grid) {
          z_minor_grid_ = z_minor_grid;
      }
  
      bool axes_type::axes_aspect_ratio_auto() const {
          return axes_aspect_ratio_auto_;
      }
  
      void axes_type::axes_aspect_ratio_auto(bool axes_aspect_ratio_auto) {
          axes_aspect_ratio_auto_ = axes_aspect_ratio_auto;
          touch();
      }
  
      void axes_type::axes_aspect_ratio(float axes_aspect_ratio) {
          axes_aspect_ratio_ = axes_aspect_ratio;
          axes_aspect_ratio_auto_ = false;
          touch();
      }
  
      float axes_type::axes_aspect_ratio() const { return axes_aspect_ratio_; }
  
      float axes_type::azimuth() const { return azimuth_; }
  
      void axes_type::azimuth(float azimuth) {
          if (azimuth != azimuth_) {
              azimuth_ = azimuth;
              if (parent_->children_.size() == 1 && !parent_->quiet_mode()) {
                  // instead of touch(), we just run the view command and replot
                  float gnuplot_azimuth_ = azimuth_ + 360;
                  while (gnuplot_azimuth_ < 0) {
                      gnuplot_azimuth_ += 360;
                  }
                  while (gnuplot_azimuth_ > 360) {
                      gnuplot_azimuth_ -= 360;
                  }
                  run_command("set view ," + num2str(gnuplot_azimuth_));
                  run_command("replot");
                  parent_->flush_commands();
              } else {
                  touch();
              }
          }
      }
  
      float axes_type::elevation() const { return elevation_; }
  
      void axes_type::elevation(float elevation) {
          if (elevation != elevation_) {
              elevation_ = elevation;
              if (parent_->children_.size() == 1 && !parent_->quiet_mode()) {
                  // instead of touch(), we just run the view command and replot
                  float gnuplot_elevation_ = 90 - elevation_;
                  while (gnuplot_elevation_ < 0) {
                      gnuplot_elevation_ += 180;
                  }
                  while (gnuplot_elevation_ > 180) {
                      gnuplot_elevation_ -= 180;
                  }
                  run_command("set view " + num2str(gnuplot_elevation_));
                  run_command("replot");
                  parent_->flush_commands();
              } else {
                  touch();
              }
          }
      }
  
      std::pair<float, float> axes_type::view() const {
          return std::make_pair(azimuth_, elevation_);
      }
  
      void axes_type::view(float azimuth, float elevation) {
          if (azimuth != azimuth_ || elevation != elevation_) {
              azimuth_ = azimuth;
              elevation_ = elevation;
              if (parent_->children_.size() == 1 && !parent_->quiet_mode()) {
                  // instead of touch(), we just run the view command and replot
                  float gnuplot_azimuth_ = azimuth_ + 360;
                  while (gnuplot_azimuth_ < 0) {
                      gnuplot_azimuth_ += 360;
                  }
                  while (gnuplot_azimuth_ > 360) {
                      gnuplot_azimuth_ -= 360;
                  }
  
                  float gnuplot_elevation_ = 90 - elevation_;
                  while (gnuplot_elevation_ < 0) {
                      gnuplot_elevation_ += 180;
                  }
                  while (gnuplot_elevation_ > 180) {
                      gnuplot_elevation_ -= 180;
                  }
  
                  run_command("set view " + num2str(gnuplot_elevation_) + "," +
                              num2str(gnuplot_azimuth_));
                  const bool is_map = azimuth_ == 0 && elevation_ == 90;
                  if (is_map) {
                      touch();
                  } else {
                      run_command("replot");
                      parent_->flush_commands();
                  }
              } else {
                  touch();
              }
          }
      }
  
      void axes_type::rotate(float azimuth, float elevation) {
          azimuth_ += azimuth;
          elevation_ += elevation;
          touch();
      }
  
      const std::vector<std::vector<double>> &axes_type::colormap() const {
          return colormap_;
      }
  
      color_array axes_type::colormap_interpolation(double value, double min,
                                                    double max) {
          return ::matplot::colormap_interpolation(value, min, max, colormap_);
      }
  
      void axes_type::colormap(const std::vector<std::vector<double>> &colormap) {
          colormap_ = colormap;
          touch();
      }
  
      bool axes_type::color_box() const { return cb_axis_.visible(); }
  
      void axes_type::color_box(bool v) {
          cb_axis_.visible(v);
          touch();
      }
  
      bool axes_type::color_box_log_scale() const {
          return cb_axis_.scale() == axis_type::axis_scale::log;
      }
  
      void axes_type::color_box_log_scale(bool v) {
          cb_axis_.scale(v ? axis_type::axis_scale::log : axis_type::axis_scale::linear);
          touch();
      }
  
      bool axes_type::grid_front() const { return grid_front_; }
  
      void axes_type::grid_front(bool grid_front) {
          grid_front_ = grid_front;
          touch();
      }
  
      bool axes_type::r_grid() const { return r_grid_; }
  
      void axes_type::r_grid(bool r_grid) {
          r_grid_ = r_grid;
          r_user_grid_ = true;
          touch();
      }
  
      bool axes_type::r_minor_grid() const { return r_minor_grid_; }
  
      void axes_type::r_minor_grid(bool r_minor_grid) {
          r_minor_grid_ = r_minor_grid;
          r_user_grid_ = true;
          touch();
      }
  
      const class axis_type &axes_type::t_axis() const { return t_axis_; }
  
      class axis_type &axes_type::t_axis() {
          return t_axis_;
      }
  
      void axes_type::t_axis(const class axis_type &t_axis) {
          t_axis_ = t_axis;
          touch();
      }
  
      size_t axes_type::max_colors() const { return max_colors_; }
  
      void axes_type::max_colors(size_t max_colors) {
          max_colors_ = max_colors;
          touch();
      }
  
      const std::pair<double, double> axes_type::color_box_range() const {
          auto a = cb_axis_.limits();
          return std::make_pair(a[0], a[1]);
      }
  
      void axes_type::color_box_range(const std::pair<double, double> &r) {
          cb_axis_.limits({r.first, r.second});
          touch();
      }
  
      void axes_type::color_box_range(double cb_min, double cb_max) {
          color_box_range({cb_min, cb_max});
      }
  
      const std::vector<std::string> &axes_type::xticklabels() const {
          return x_axis().ticklabels();
      }
  
      void axes_type::xticklabels(const std::vector<std::string> &ticks) {
          x_axis().ticklabels(ticks);
          touch();
      }
  
      const std::vector<std::string> &axes_type::x2ticklabels() const {
          return x2_axis().ticklabels();
      }
  
      void axes_type::x2ticklabels(const std::vector<std::string> &ticks) {
          x2_axis().ticklabels(ticks);
          touch();
      }
  
      const std::vector<std::string> &axes_type::yticklabels() const {
          return y_axis().ticklabels();
      }
  
      void axes_type::yticklabels(const std::vector<std::string> &ticks) {
          y_axis().ticklabels(ticks);
          touch();
      }
  
      const std::vector<std::string> &axes_type::y2ticklabels() const {
          return y2_axis().ticklabels();
      }
  
      void axes_type::y2ticklabels(const std::vector<std::string> &ticks) {
          y2_axis().ticklabels(ticks);
          touch();
      }
  
      const std::vector<std::string> &axes_type::zticklabels() const {
          return z_axis().ticklabels();
      }
  
      void axes_type::zticklabels(const std::vector<std::string> &ticks) {
          z_axis().ticklabels(ticks);
          touch();
      }
  
      const std::vector<std::string> &axes_type::rticklabels() const {
          return r_axis().ticklabels();
      }
  
      void axes_type::rticklabels(const std::vector<std::string> &ticks) {
          r_axis().ticklabels(ticks);
          touch();
      }
  
      void axes_type::xtickangle(double degrees) {
          x_axis().tickangle(static_cast<float>(degrees));
      }
  
      void axes_type::x2tickangle(double degrees) {
          x2_axis().tickangle(static_cast<float>(degrees));
      }
  
      void axes_type::ytickangle(double degrees) {
          y_axis().tickangle(static_cast<float>(degrees));
      }
  
      void axes_type::y2tickangle(double degrees) {
          y2_axis().tickangle(static_cast<float>(degrees));
      }
  
      void axes_type::ztickangle(double degrees) {
          z_axis().tickangle(static_cast<float>(degrees));
      }
  
      void axes_type::rtickangle(double degrees) {
          r_axis().tickangle(static_cast<float>(degrees));
      }
  
      double axes_type::xtickangle() { return x_axis().tickangle(); }
  
      double axes_type::x2tickangle() { return x2_axis().tickangle(); }
  
      double axes_type::ytickangle() { return y_axis().tickangle(); }
  
      double axes_type::y2tickangle() { return y2_axis().tickangle(); }
  
      double axes_type::ztickangle() { return z_axis().tickangle(); }
  
      double axes_type::rtickangle() { return r_axis().tickangle(); }
  
      const color_array &axes_type::title_color() const { return title_color_; }
  
      void axes_type::title_color(const color_array &title_color) {
          title_color_ = title_color;
          touch();
      }
  
      bool axes_type::title_enhanced() const { return title_enhanced_; }
  
      void axes_type::title_enhanced(bool title_enhanced) {
          title_enhanced_ = title_enhanced;
          touch();
      }
  
      std::array<double, 2> axes_type::xlim() const {
          if (x_axis_.limits_mode_auto()) {
              auto [xmin, xmax, ymin, ymax] = this->child_limits();
              return std::array<double, 2>{xmin, xmax};
          } else {
              return x_axis().limits();
          }
      }
  
      void axes_type::xlim(const std::array<double, 2> &lim) {
          x_axis().limits(lim);
          touch();
      }
  
      std::array<double, 2> axes_type::x2lim() const {
          return x2_axis().limits();
      }
  
      void axes_type::x2lim(const std::array<double, 2> &lim) {
          x2_axis().limits(lim);
          touch();
      }
  
      std::array<double, 2> axes_type::ylim() const {
          if (x_axis_.limits_mode_auto()) {
              auto [xmin, xmax, ymin, ymax] = this->child_limits();
              return std::array<double, 2>{ymin, ymax};
          } else {
              return y_axis().limits();
          }
      }
  
      void axes_type::ylim(const std::array<double, 2> &lim) {
          y_axis().limits(lim);
          touch();
      }
  
      std::array<double, 2> axes_type::y2lim() const {
          return y2_axis().limits();
      }
  
      void axes_type::y2lim(const std::array<double, 2> &lim) {
          y2_axis().limits(lim);
          touch();
      }
  
      std::array<double, 2> axes_type::zlim() const { return z_axis().limits(); }
  
      void axes_type::zlim(const std::array<double, 2> &lim) {
          z_axis().limits(lim);
          touch();
      }
  
      std::array<double, 2> axes_type::rlim() const { return r_axis().limits(); }
  
      void axes_type::rlim(const std::array<double, 2> &lim) {
          r_axis().limits(lim);
          touch();
      }
  
      std::array<double, 2> axes_type::tlim() const { return t_axis().limits(); }
  
      void axes_type::tlim(const std::array<double, 2> &lim) {
          t_axis().limits(lim);
          touch();
      }
  
      void axes_type::xlim(keyword_automatic_type) {
          x_axis().limits_mode_manual(false);
      }
  
      void axes_type::xlim(keyword_manual_type) {
          x_axis().limits_mode_manual(true);
      }
  
      void axes_type::x2lim(keyword_automatic_type) {
          x2_axis().limits_mode_manual(false);
      }
  
      void axes_type::x2lim(keyword_manual_type) {
          x2_axis().limits_mode_manual(true);
      }
  
      void axes_type::ylim(keyword_automatic_type) {
          y_axis().limits_mode_manual(false);
      }
  
      void axes_type::ylim(keyword_manual_type) {
          y_axis().limits_mode_manual(true);
      }
  
      void axes_type::y2lim(keyword_automatic_type) {
          y2_axis().limits_mode_manual(false);
      }
  
      void axes_type::y2lim(keyword_manual_type) {
          y2_axis().limits_mode_manual(true);
      }
  
      void axes_type::zlim(keyword_automatic_type) {
          z_axis().limits_mode_manual(false);
      }
  
      void axes_type::zlim(keyword_manual_type) {
          z_axis().limits_mode_manual(true);
      }
  
      void axes_type::rlim(keyword_automatic_type) {
          r_axis().limits_mode_manual(false);
      }
  
      void axes_type::rlim(keyword_manual_type) {
          r_axis().limits_mode_manual(true);
      }
  
      void axes_type::tlim(keyword_automatic_type) {
          t_axis().limits_mode_manual(false);
      }
  
      void axes_type::tlim(keyword_manual_type) {
          t_axis().limits_mode_manual(true);
      }
  
      void axes_type::clear() {
          children_.clear();
          legend_.reset();
          colororder_index_ = 0;
          x_axis_.limits_mode_auto(true);
          y_axis_.limits_mode_auto(true);
          y2_axis_.limits_mode_auto(true);
          z_axis_.limits_mode_auto(true);
          touch();
      }
  
      const class axis_type &axes_type::cb_axis() const { return cb_axis_; }
  
      class axis_type &axes_type::cb_axis() {
          return cb_axis_;
      }
  
      std::array<double, 2> axes_type::cblim() const { return cb_axis_.limits(); }
  
      void axes_type::cblim(const std::array<double, 2> &a) {
          cb_axis_.limits(a);
      }
  
      void axes_type::cblim(keyword_automatic_type) {
          cb_axis_.limits_mode_auto(true);
      }
  
      void axes_type::cblim(keyword_manual_type) {
          cb_axis_.limits_mode_manual(true);
      }
  
      const std::string &axes_type::cblabel() const { return cb_axis_.label(); }
  
      void axes_type::cblabel(std::string_view str) { cb_axis_.label(); }
  
      const std::string &axes_type::cbtickformat() const {
          return cb_axis_.tick_label_format();
      }
  
      void axes_type::cbtickformat(std::string_view str) {
          cb_axis_.tick_label_format(str);
      }
  
      const std::vector<double> &axes_type::cbticks() const {
          return cb_axis_.tick_values();
      }
  
      void axes_type::cbticks(const std::vector<double> &ticks) {
          cb_axis_.tick_values(ticks);
      }
  
      const std::vector<std::string> &axes_type::cbticklabels() const {
          return cb_axis_.ticklabels();
      }
  
      void axes_type::cbticklabels(const std::vector<std::string> &labels) {
          cb_axis_.ticklabels();
      }
  
      void axes_type::cbtickangle(double degrees) { cb_axis_.tickangle(); }
  
      double axes_type::cbtickangle() { return cb_axis_.tickangle(); }
  
      bool axes_type::cb_vertical() const { return cb_vertical_; }
  
      void axes_type::cb_vertical(bool cb_vertical) {
          cb_vertical_ = cb_vertical;
          touch();
      }
  
      const std::array<float, 4> &axes_type::cb_position() const {
          return cb_position_;
      }
  
      void axes_type::cb_position(const std::array<float, 4> &cb_position) {
          cb_position_ = cb_position;
          cb_vertical_ = cb_position_[2] <= cb_position_[3];
          touch();
      }
  
      bool axes_type::cb_inside() const { return cb_inside_; }
  
      void axes_type::cb_inside(bool cb_inside) {
          cb_inside_ = cb_inside;
          touch();
      }
  
      class axes_silencer {
        public:
          explicit axes_silencer(class axes_type *ax)
              : ax_(ax), previous_quiet_mode_(ax->parent()->quiet_mode()) {
              ax_->parent()->quiet_mode(true);
          }
  
          ~axes_silencer() {
              ax_->parent()->quiet_mode(previous_quiet_mode_);
              if (!previous_quiet_mode_) {
                  ax_->draw();
              }
          }
  
        private:
          class axes_type *ax_;
  
          bool previous_quiet_mode_;
      };
  
      line_handle axes_type::plot(const std::vector<double> &x,
                                  const std::vector<double> &y,
                                  std::string_view line_spec) {
          axes_silencer s{this};
          line_handle l = std::make_shared<class line>(this, x, y, line_spec);
          this->emplace_object(l);
          return l;
      }
  
      line_handle axes_type::plot(const std::vector<double> &y,
                                  std::string_view line_spec) {
          axes_silencer s{this};
          line_handle l = std::make_shared<class line>(this, y, line_spec);
          this->emplace_object(l);
          return l;
      }
  
      std::vector<line_handle>
      axes_type::plot(const std::vector<double> &x,
                      const std::vector<std::vector<double>> &Y,
                      std::string_view line_spec) {
          axes_silencer s{this};
          std::vector<line_handle> hs;
          bool previous_replace_status = this->next_plot_replace();
          for (const auto &y : Y) {
              hs.emplace_back(this->plot(x, y, line_spec));
              this->next_plot_replace(false);
          }
          this->next_plot_replace(previous_replace_status);
          return hs;
      }
  
      std::vector<line_handle>
      axes_type::plot(const std::vector<std::vector<double>> &Y,
                      std::string_view line_spec) {
          axes_silencer s{this};
          std::vector<line_handle> hs;
          bool previous_replace_status = this->next_plot_replace();
          for (const auto &y : Y) {
              hs.emplace_back(this->plot(y, line_spec));
              this->next_plot_replace(false);
          }
          this->next_plot_replace(previous_replace_status);
          return hs;
      }
  
      std::vector<line_handle>
      axes_type::rgbplot(const std::vector<std::vector<double>> &map) {
          axes_silencer s{this};
          // Process data
          auto rgb = transpose(map);
          auto r = plot(rgb[0], "r-");
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
          auto g = plot(rgb[1], "g-");
          auto b = plot(rgb[2], "b-");
          // If nothing uses the palette, we can't show it
          // Create markers using it on the first points
          double n = static_cast<double>(map.size());
          auto palette_enabler =
              plot(vector_1d{1, 1, 1, n, n, n},
                   vector_1d{rgb[0].front(), rgb[1].front(), rgb[2].front(),
                             rgb[0].back(), rgb[1].back(), rgb[2].back()},
                   ".");
          palette_enabler->marker_colors({0., 0., 0., 1., 1., 1.});
          this->next_plot_replace(p2);
          return {r, g, b, palette_enabler};
      }
  
      line_handle axes_type::plot3(const std::vector<double> &x,
                                   const std::vector<double> &y,
                                   const std::vector<double> &z,
                                   std::string_view line_spec) {
          axes_silencer s{this};
          line_handle l = std::make_shared<class line>(this, x, y, z, line_spec);
          this->emplace_object(l);
          return l;
      }
  
      std::vector<line_handle>
      axes_type::plot3(const std::vector<std::vector<double>> &X,
                       const std::vector<std::vector<double>> &Y,
                       const std::vector<double> &z,
                       std::string_view line_spec) {
          axes_silencer s{this};
          auto it_x = X.begin();
          auto it_y = Y.begin();
          bool previous_status = this->next_plot_replace();
          this->next_plot_replace(false);
          std::vector<line_handle> ls;
          while (it_x != X.end() && it_y != Y.end()) {
              ls.emplace_back(this->plot3(*it_x, *it_y, z, line_spec));
              ++it_x;
              ++it_y;
          }
          this->next_plot_replace(previous_status);
          return ls;
      }
  
      std::vector<line_handle>
      axes_type::plot3(const std::vector<std::vector<double>> &X,
                       const std::vector<std::vector<double>> &Y,
                       const std::vector<std::vector<double>> &Z,
                       std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          auto it_x = X.begin();
          auto it_y = Y.begin();
          auto it_z = Z.begin();
          bool previous_status = this->next_plot_replace();
          this->next_plot_replace(false);
          std::vector<line_handle> ls;
          while (it_x != X.end() && it_y != Y.end() && it_z != Z.end()) {
              ls.emplace_back(this->plot3(*it_x, *it_y, *it_z, line_spec));
              ++it_x;
              ++it_y;
              ++it_z;
          }
          this->next_plot_replace(previous_status);
          return ls;
      }
  
      stair_handle axes_type::stairs(const std::vector<double> &x,
                                     const std::vector<double> &y,
                                     std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          stair_handle l = std::make_shared<class stair>(this, x, y, line_spec);
          this->emplace_object(l);
          return l;
      }
  
      stair_handle axes_type::stairs(const std::vector<double> &y,
                                     std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          stair_handle l = std::make_shared<class stair>(this, y, line_spec);
          this->emplace_object(l);
          return l;
      }
  
      std::vector<stair_handle>
      axes_type::stairs(const std::vector<std::vector<double>> &X,
                        const std::vector<std::vector<double>> &Y,
                        std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          std::vector<stair_handle> hs;
          bool previous_replace_status = this->next_plot_replace();
          auto X_it = X.begin();
          for (const auto &y : Y) {
              hs.emplace_back(this->stairs(*X_it, y, line_spec));
              ++X_it;
              this->next_plot_replace(false);
          }
          this->next_plot_replace(previous_replace_status);
          return hs;
      }
  
      std::vector<stair_handle>
      axes_type::stairs(const std::vector<std::vector<double>> &Y,
                        std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          std::vector<stair_handle> hs;
          bool previous_replace_status = this->next_plot_replace();
          for (const auto &y : Y) {
              hs.emplace_back(this->stairs(y, line_spec));
              this->next_plot_replace(false);
          }
          this->next_plot_replace(previous_replace_status);
          return hs;
      }
  
      std::vector<stair_handle>
      axes_type::stairs(const std::vector<double> &x,
                        const std::vector<std::vector<double>> &Y,
                        std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          std::vector<stair_handle> hs;
          bool previous_replace_status = this->next_plot_replace();
          for (const auto &y : Y) {
              hs.emplace_back(this->stairs(x, y, line_spec));
              this->next_plot_replace(false);
          }
          this->next_plot_replace(previous_replace_status);
          return hs;
      }
  
      error_bar_handle axes_type::errorbar(const std::vector<double> &x,
                                           const std::vector<double> &y,
                                           const std::vector<double> &error,
                                           error_bar::type type,
                                           std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          error_bar_handle l = std::make_shared<class error_bar>(
              this, x, y, error, type, line_spec);
          this->emplace_object(l);
          return l;
      }
  
      error_bar_handle axes_type::errorbar(const std::vector<double> &x,
                                           const std::vector<double> &y,
                                           const std::vector<double> &y_neg_delta,
                                           const std::vector<double> &y_pos_delta,
                                           const std::vector<double> &x_neg_delta,
                                           const std::vector<double> &x_pos_delta,
                                           std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          error_bar_handle l = std::make_shared<class error_bar>(
              this, x, y, y_neg_delta, y_pos_delta, x_neg_delta, x_pos_delta,
              line_spec);
          this->emplace_object(l);
          return l;
      }
  
      error_bar_handle axes_type::errorbar(const std::vector<double> &x,
                                           const std::vector<double> &y,
                                           const std::vector<double> &y_error,
                                           std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          return this->errorbar(x, y, y_error, error_bar::type::vertical,
                                line_spec);
      }
  
      std::vector<filled_area_handle> axes_type::area(
          const std::vector<double> &x, const std::vector<std::vector<double>> &Y,
          double base_value, bool stacked, std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          std::vector<filled_area_handle> hs;
          bool previous_state = this->next_plot_replace();
          // we are going to add y_i in reverse, so we also get the colors
          // in reverse, as if we were getting the colors in order
          std::vector<color_array> colors;
          for (size_t i = 0; i < Y.size(); ++i) {
              colors.emplace_back(this->get_color_and_bump());
              if (!stacked) {
                  // reduce opacity if not stacked
                  colors.back()[0] = 1.f - (1.f - colors.back()[0]) * 0.7f;
              }
          }
          // use reverse iterator in this case because the first Ys should be
          // the ones in the background. The reverse happens for other object
          // types
          for (auto [Y_it, c_it] = std::make_pair(Y.rbegin(), colors.rbegin());
               Y_it != Y.rend(); ++Y_it, ++c_it) {
              const auto &y = *Y_it;
              auto h = std::make_shared<class filled_area>(
                  this, x, y, std::vector({base_value}), stacked, line_spec);
              h->line_width(1.0);
              h->face_color(*c_it);
              this->emplace_object(h);
              hs.emplace_back(h);
              this->next_plot_replace(false);
          }
          this->next_plot_replace(previous_state);
          return hs;
      }
  
      filled_area_handle axes_type::area(const std::vector<double> &x,
                                         const std::vector<double> &y,
                                         double base_value, bool stacked,
                                         std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          auto as =
              this->area(x, std::vector({y}), base_value, stacked, line_spec);
          return as[0];
      }
  
      std::vector<filled_area_handle>
      axes_type::area(const std::vector<std::vector<double>> &Y,
                      double base_value, bool stacked,
                      std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          return this->area(iota(1., static_cast<double>(Y.begin()->size())), Y, base_value, stacked,
                            line_spec);
      }
  
      std::vector<filled_area_handle>
      axes_type::area(const std::vector<std::vector<double>> &Y, bool stacked,
                      std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          return this->area(iota(1., static_cast<double>(Y.begin()->size())), Y, 0., stacked,
                            line_spec);
      }
  
      filled_area_handle axes_type::area(const std::vector<double> &y,
                                         double base_value, bool stacked,
                                         std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          return this->area(iota(1., static_cast<double>(y.size())), {y},
                            base_value, stacked, line_spec);
      }
  
      line_handle axes_type::fimplicit(axes_type::implicit_function_type equation,
                                       const std::array<double, 4> &xy_interval,
                                       std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          auto [X, Y] = meshgrid(linspace(xy_interval[0], xy_interval[1]),
                                 linspace(xy_interval[2], xy_interval[3]));
          auto Z = transform(X, Y, equation);
          std::vector<std::pair<vector_1d, vector_1d>> c =
              contourc(X, Y, Z, vector_1d{0.});
          if (!c.empty()) {
              auto l = this->plot(c[0].first, c[0].second, line_spec);
              return l;
          } else {
              return this->plot(vector_1d{0.}, vector_1d{0.}, line_spec);
          }
      }
  
      line_handle axes_type::fimplicit(axes_type::implicit_function_type equation,
                                       std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          return this->fimplicit(equation, {-5, 5, -5, 5}, line_spec);
      }
  
      string_function_handle axes_type::fplot(std::string_view equation,
                                              std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          string_function_handle obj =
              std::make_shared<class string_function>(this, equation, line_spec);
          this->emplace_object(obj);
          obj->touch();
          return obj;
      }
  
      std::vector<string_function_handle>
      axes_type::fplot(std::vector<std::string> equations,
                       std::vector<std::string> line_specs) {
          axes_silencer temp_silencer_{this};
          std::vector<string_function_handle> res;
          auto it_line_specs = line_specs.begin();
          for (const auto &equation : equations) {
              if (it_line_specs != line_specs.end()) {
                  res.emplace_back(this->fplot(equation, *it_line_specs));
              } else {
                  res.emplace_back(this->fplot(equation));
              }
              ++it_line_specs;
          }
          return res;
      }
  
      function_line_handle axes_type::fplot(function_line::function_type equation,
                                            const std::array<double, 2> &x_range,
                                            std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          function_line_handle obj = std::make_shared<class function_line>(
              this, equation, x_range, line_spec);
          this->emplace_object(obj);
          obj->touch();
          return obj;
      }
  
      function_line_handle axes_type::fplot(function_line::function_type equation,
                                            std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          function_line_handle obj = std::make_shared<class function_line>(
              this, equation, std::array<double, 2>({-5, 5}), line_spec);
          this->emplace_object(obj);
          obj->touch();
          return obj;
      }
  
      function_line_handle
      axes_type::fplot(function_line::function_type function_x,
                       function_line::function_type function_y,
                       const std::array<double, 2> &t_range,
                       std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          function_line_handle obj = std::make_shared<class function_line>(
              this, function_x, function_y, t_range, line_spec);
          this->emplace_object(obj);
          obj->touch();
          return obj;
      }
  
      std::vector<function_line_handle>
      axes_type::fplot(std::vector<function_line::function_type> equations,
                       std::array<double, 2> x_range,
                       std::vector<std::string> line_specs) {
          axes_silencer temp_silencer_{this};
          std::vector<function_line_handle> res;
          auto it_line_specs = line_specs.begin();
          for (const auto &equation : equations) {
              if (it_line_specs != line_specs.end()) {
                  res.emplace_back(
                      this->fplot(equation, x_range, *it_line_specs));
              } else {
                  res.emplace_back(this->fplot(equation, x_range));
              }
              ++it_line_specs;
          }
          return res;
      }
  
      std::vector<function_line_handle>
      axes_type::fplot(std::vector<function_line::function_type> equations,
                       std::vector<double> x_range,
                       std::vector<std::string> line_specs) {
          return this->fplot(equations, to_array<2>(x_range), line_specs);
      }
  
      function_line_handle
      axes_type::fplot3(function_line::function_type function_x,
                        function_line::function_type function_y,
                        function_line::function_type function_z,
                        const std::array<double, 2> &t_range,
                        std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          function_line_handle obj = std::make_shared<class function_line>(
              this, function_x, function_y, function_z, t_range, line_spec);
          this->emplace_object(obj);
          obj->touch();
          return obj;
      }
  
      /// Histogram - Choose binning algorithm and normalization algorithm
      histogram_handle
      axes_type::hist(const std::vector<double> &data,
                      histogram::binning_algorithm algorithm,
                      enum histogram::normalization normalization_alg) {
          axes_silencer temp_silencer_{this};
          histogram_handle l = std::make_shared<class histogram>(
              this, data, algorithm, normalization_alg);
          this->emplace_object(l);
          return l;
      }
  
      /// Histogram - Fixed number of bins
      /// If n_bins = 0, automatic number of bins
      histogram_handle axes_type::hist(const std::vector<double> &data,
                                       size_t n_bins) {
          axes_silencer temp_silencer_{this};
          histogram_handle l =
              std::make_shared<class histogram>(this, data, n_bins);
          this->emplace_object(l);
          return l;
      }
  
      /// Histogram - Fixed edges
      histogram_handle axes_type::hist(const std::vector<double> &data,
                                       const std::vector<double> &edges) {
          axes_silencer temp_silencer_{this};
          histogram_handle l =
              std::make_shared<class histogram>(this, data, edges);
          this->emplace_object(l);
          return l;
      }
  
      /// Histogram - Normalization algorithm
      histogram_handle
      axes_type::hist(const std::vector<double> &data,
                      enum histogram::normalization normalization_alg) {
          axes_silencer temp_silencer_{this};
          return this->hist(data, histogram::binning_algorithm::automatic,
                            normalization_alg);
      }
  
      /// Histogram - Categorical histogram from strings
      histogram_handle
      axes_type::hist(const std::vector<std::string> &data,
                      enum histogram::normalization normalization_alg) {
          axes_silencer temp_silencer_{this};
          std::map<std::string, double> category_indexes;
          std::vector<double> numeric_data;
          for (const auto &category : data) {
              auto it = category_indexes.find(category);
              if (it != category_indexes.end()) {
                  numeric_data.emplace_back(it->second);
              } else {
                  double category_num = static_cast<double>(category_indexes.size()) + 1;
                  category_indexes[category] = category_num;
                  numeric_data.emplace_back(category_num);
              }
          }
          std::vector<double> fixed_edges =
              iota(0.5, 1, 0.5 + category_indexes.size());
          histogram_handle l = std::make_shared<class histogram>(
              this, numeric_data, fixed_edges, normalization_alg);
          this->emplace_object(l);
          this->xticks(iota(1., static_cast<double>(category_indexes.size())));
  
          std::vector<std::string> unique_categories(category_indexes.size());
          for (const auto &[category, index] : category_indexes) {
              unique_categories[static_cast<size_t>(round(index) - 1)] = category;
          }
          this->x_axis().ticklabels(unique_categories);
          l->bar_width(0.9);
          this->x_axis().limits(
              {0.5, static_cast<double>(category_indexes.size()) + 0.5});
          return l;
      }
  
      axes_object_handle axes_type::binscatter(
          const std::vector<double> &x, const std::vector<double> &y,
          const std::vector<double> &edges_x, const std::vector<double> &edges_y,
          enum bin_scatter_style scatter_style,
          enum histogram::normalization normalization_alg) {
          axes_silencer temp_silencer_{this};
  
          // count data in each bin combination
          std::vector<std::vector<double>> values =
              histcounts2(x, y, edges_x, edges_y, normalization_alg);
          double max_value = values[0][0];
          for (const auto &row : values) {
              for (const auto &value : row) {
                  if (value > max_value) {
                      max_value = value;
                  }
              }
          }
  
          // set automatic mode
          if (scatter_style == bin_scatter_style::automatic) {
              scatter_style = bin_scatter_style::point_size;
          }
  
          // value stats
          double sum_values = normalization_alg == histogram::normalization::count
                                  ? x.size()
                                  : 0.;
          if (normalization_alg != histogram::normalization::count) {
              for (const auto &row : values) {
                  for (const auto &v : row) {
                      sum_values += v;
                  }
              }
          }
  
          // stats for jitter style
          double n_squares = static_cast<double>(values.size() * values[0].size());
          double avg_square_value = sum_values / n_squares;
          double n_point_per_plot = 500.;
          double avg_pts_per_square = n_point_per_plot / n_squares;
  
          // the point sizes must be proportional to the max value and the
          // area available for each point
          double x_bin_width =
              edges_x[edges_x.size() / 2 + 1] - edges_x[edges_x.size() / 2];
          double y_bin_width =
              edges_y[edges_y.size() / 2 + 1] - edges_y[edges_y.size() / 2];
          double reference_width = std::max(x_bin_width, y_bin_width);
  
          // create bin positions and their sizes
          std::vector<double> bin_x;
          std::vector<double> bin_y;
          std::vector<double> point_sizes;
          std::vector<double> colors;
          for (size_t i = 0; i < edges_x.size() - 1; ++i) {
              for (size_t j = 0; j < edges_y.size() - 1; ++j) {
                  if (values[i][j] != 0.0) {
                      if (scatter_style != bin_scatter_style::jitter) {
                          bin_x.emplace_back((edges_x[i] + edges_x[i + 1]) / 2.);
                          bin_y.emplace_back((edges_y[j] + edges_y[j + 1]) / 2.);
                      }
                      switch (scatter_style) {
                      case bin_scatter_style::point_size:
                          point_sizes.emplace_back((values[i][j] / max_value) *
                                                   reference_width * 20);
                          break;
                      case bin_scatter_style::point_alpha:
                      case bin_scatter_style::point_colormap:
                      case bin_scatter_style::heatmap:
                          colors.emplace_back(values[i][j]);
                          break;
                      case bin_scatter_style::jitter: {
                          double n_expected_points_in_bin = values[i][j] /
                                                            avg_square_value *
                                                            avg_pts_per_square;
                          double integer_part = floor(n_expected_points_in_bin);
                          double fractional_part =
                              n_expected_points_in_bin - integer_part;
                          size_t n_points_in_bin =
                              static_cast<size_t>(integer_part + (rand(0, 1) < fractional_part));
                          for (size_t k = 0; k < n_points_in_bin; ++k) {
                              bin_x.emplace_back(
                                  rand(edges_x[i], edges_x[i + 1]));
                              bin_y.emplace_back(
                                  rand(edges_y[j], edges_y[j + 1]));
                          }
                          break;
                      }
                      case bin_scatter_style::automatic:
                          break;
                      }
                  }
              }
          }
          // create scatter plot with this data
          axes_object_handle r;
          line_handle s;
          matrix_handle m;
          switch (scatter_style) {
          case bin_scatter_style::point_size: {
              s = this->scatter(bin_x, bin_y, point_sizes);
              auto c = s->marker_color();
              c[0] = 0.2;
              s->marker_color(c);
              s->marker_face(true);
              r = s;
              break;
          }
          case bin_scatter_style::point_alpha: {
              s = this->scatter(bin_x, bin_y, std::vector<double>{}, colors);
              auto c = s->marker_color();
              std::vector<double> alpha_100 = {1, 1, 1};
              std::vector<double> alpha_0 = {c[1] - c[1] / 3., c[2] - c[2] / 3.,
                                             c[3] - c[3] / 3.};
              this->colormap({alpha_100, alpha_0});
              s->marker_face(true);
              this->color_box(true);
              this->color_box_log_scale(true);
              r = s;
              break;
          }
          case bin_scatter_style::heatmap:
          case bin_scatter_style::automatic:
              this->colormap(palette::blues());
              m = std::make_shared<class matrix>(this, values);
              m->always_hide_labels(true);
              m->x(edges_x.front());
              m->w(edges_x.back() - edges_x.front());
              m->y(edges_y.front());
              m->h(edges_y.back() - edges_y.front());
              this->emplace_object(m);
              this->color_box(true);
              this->y_axis().reverse(false);
              this->colormap(palette::blues());
              r = m;
              break;
          case bin_scatter_style::point_colormap:
              s = this->scatter(bin_x, bin_y, std::vector<double>{}, colors);
              s->marker_face(true);
              this->color_box(true);
              this->color_box_log_scale(true);
              r = s;
              break;
          case bin_scatter_style::jitter:
              s = this->scatter(bin_x, bin_y);
              r = s;
              break;
          }
  
          return r;
      }
  
      axes_object_handle
      axes_type::binscatter(const std::vector<double> &x,
                            const std::vector<double> &y,
                            histogram::binning_algorithm algorithm,
                            enum bin_scatter_style scatter_style,
                            enum histogram::normalization normalization_alg) {
          auto [min_x, max_x] = std::minmax_element(x.begin(), x.end());
          auto x_edges =
              histogram::histogram_edges(x, *min_x, *max_x, algorithm, false);
  
          auto [min_y, max_y] = std::minmax_element(y.begin(), y.end());
          auto y_edges =
              histogram::histogram_edges(y, *min_y, *max_y, algorithm, false);
  
          return this->binscatter(x, y, x_edges, y_edges, scatter_style,
                                  normalization_alg);
      }
  
      axes_object_handle
      axes_type::binscatter(const std::vector<double> &x,
                            const std::vector<double> &y,
                            enum bin_scatter_style scatter_style,
                            histogram::binning_algorithm algorithm,
                            enum histogram::normalization normalization_alg) {
          return this->binscatter(x, y, algorithm, scatter_style,
                                  normalization_alg);
      }
  
      axes_object_handle
      axes_type::binscatter(const std::vector<double> &x,
                            const std::vector<double> &y, size_t nbins_x,
                            size_t nbins_y, enum bin_scatter_style scatter_style,
                            enum histogram::normalization normalization_alg) {
          auto [minx, maxx] = std::minmax_element(x.begin(), x.end());
          double xrange = *maxx - *minx;
          auto x_edges =
              histogram::bin_picker(*minx, *maxx, nbins_x, xrange / nbins_x);
  
          auto [miny, maxy] = std::minmax_element(y.begin(), y.end());
          double yrange = *maxy - *miny;
          auto y_edges =
              histogram::bin_picker(*miny, *maxy, nbins_y, yrange / nbins_y);
  
          return this->binscatter(x, y, x_edges, y_edges, scatter_style,
                                  normalization_alg);
      }
  
      bars_handle axes_type::bar(const std::vector<double> &x,
                                 const std::vector<double> &y) {
          axes_silencer temp_silencer_{this};
          bars_handle l = std::make_shared<class bars>(this, x, y);
          this->emplace_object(l);
          return l;
      }
  
      bars_handle axes_type::bar(const std::vector<double> &y) {
          axes_silencer temp_silencer_{this};
          bars_handle l = std::make_shared<class bars>(this, y);
          this->emplace_object(l);
          return l;
      }
  
      bars_handle axes_type::bar(const std::vector<std::vector<double>> &Y) {
          axes_silencer temp_silencer_{this};
          bars_handle l = std::make_shared<class bars>(this, Y);
          this->emplace_object(l);
          return l;
      }
  
      bars_handle axes_type::bar(const std::vector<double> &x,
                                 const std::vector<std::vector<double>> &Y) {
          axes_silencer temp_silencer_{this};
          bars_handle l = std::make_shared<class bars>(this, x, Y);
          this->emplace_object(l);
          return l;
      }
  
      bars_handle axes_type::bar(const std::vector<double> &y, double width) {
          axes_silencer temp_silencer_{this};
          bars_handle l = std::make_shared<class bars>(this, y);
          l->bar_width(static_cast<float>(width));
          this->emplace_object(l);
          return l;
      }
  
      /// Bar stacked - Many Y - Automatic x
      std::vector<bars_handle>
      axes_type::barstacked(const std::vector<std::vector<double>> &Y) {
          axes_silencer temp_silencer_{this};
          std::vector<bars_handle> hs;
          auto stacked_positive_values = Y[0];
          for (size_t i = 0; i < Y.size(); ++i) {
              if (Y[0][i] < 0) {
                  stacked_positive_values[i] = 0;
              }
          }
          auto stacked_negative_values = Y[0];
          for (size_t i = 0; i < Y.size(); ++i) {
              if (Y[0][i] > 0) {
                  stacked_negative_values[i] = 0;
              }
          }
          for (size_t i = 0; i < Y.size(); ++i) {
              auto current_stack = stacked_positive_values;
              for (size_t j = 0; j < Y[i].size(); ++j) {
                  if (Y[i][j] < 0) {
                      current_stack[j] = stacked_negative_values[j];
                  }
              }
              bars_handle h = std::make_shared<class bars>(this, current_stack);
              h->face_color(this->get_color_and_bump());
              hs.emplace_back(h);
              if (i != Y.size() - 1) {
                  for (size_t j = 0; j < Y[i + 1].size(); ++j) {
                      if (Y[i + 1][j] > 0) {
                          stacked_positive_values[j] += Y[i + 1][j];
                      } else {
                          stacked_negative_values[j] += Y[i + 1][j];
                      }
                  }
              }
          }
  
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
          for (auto it = hs.rbegin(); it != hs.rend(); ++it) {
              this->emplace_object(*it);
          }
          this->next_plot_replace(p2);
          return hs;
      }
  
      /// Bar stacked - One x - Many Y
      std::vector<bars_handle>
      axes_type::barstacked(const std::vector<double> &x,
                            const std::vector<std::vector<double>> &Y) {
          axes_silencer temp_silencer_{this};
          std::vector<bars_handle> hs;
          auto stacked_positive_values = Y[0];
          for (size_t i = 0; i < Y.size(); ++i) {
              if (Y[0][i] < 0) {
                  stacked_positive_values[i] = 0;
              }
          }
          auto stacked_negative_values = Y[0];
          for (size_t i = 0; i < Y.size(); ++i) {
              if (Y[0][i] > 0) {
                  stacked_negative_values[i] = 0;
              }
          }
          for (size_t i = 0; i < Y.size(); ++i) {
              auto current_stack = stacked_positive_values;
              for (size_t j = 0; j < Y[i].size(); ++j) {
                  if (Y[i][j] < 0) {
                      current_stack[j] = stacked_negative_values[j];
                  }
              }
              bars_handle h =
                  std::make_shared<class bars>(this, x, current_stack);
              h->face_color(this->get_color_and_bump());
              hs.emplace_back(h);
              if (i != Y.size() - 1) {
                  for (size_t j = 0; j < Y[i + 1].size(); ++j) {
                      if (Y[i + 1][j] > 0) {
                          stacked_positive_values[j] += Y[i + 1][j];
                      } else {
                          stacked_negative_values[j] += Y[i + 1][j];
                      }
                  }
              }
          }
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
          for (auto it = hs.rbegin(); it != hs.rend(); ++it) {
              this->emplace_object(*it);
          }
          this->next_plot_replace(p2);
          return hs;
      }
  
      matrix_handle
      axes_type::heatmap(const std::vector<std::vector<double>> &m) {
          axes_silencer temp_silencer_{this};
          // create matrix in the xlim
          matrix_handle heatmap = std::make_shared<class matrix>(this, m);
          this->emplace_object(heatmap);
          this->y_axis().reverse(true);
          this->color_box(true);
          auto [h_, w_] = size(m);
          this->y_axis().tick_values(iota(1., static_cast<double>(h_)));
          this->y_axis().tick_length(0.);
          this->x_axis().tick_values(iota(1., static_cast<double>(w_)));
          this->x_axis().tick_length(0.);
          this->colormap(palette::blues());
          this->box(true);
          this->grid(false);
          this->minor_grid(true);
          this->grid_front(true);
          this->minor_grid_line_style(line_spec("-k"));
          return heatmap;
      }
  
      matrix_handle axes_type::pcolor(const std::vector<std::vector<double>> &m) {
          axes_silencer temp_silencer_{this};
          // create matrix in the xlim
          matrix_handle pcolor = std::make_shared<class matrix>(this, m);
          pcolor->x(1.5);
          pcolor->y(1.5);
          pcolor->always_hide_labels(true);
          this->emplace_object(pcolor);
  
          auto [h, w] = size(m);
          this->ylim({1, static_cast<double>(h + 1)});
          this->y_axis().tick_values(iota(1, h));
          this->xlim({1, static_cast<double>(w + 1)});
          this->x_axis().tick_values(iota(1, w));
          this->y_axis().reverse(false);
          this->box(true);
          this->grid(true);
          this->grid_color({0., 0., 0., 0.});
          this->minor_grid(false);
          this->grid_front(true);
          this->minor_grid_line_style(line_spec("-k"));
          this->color_box(false);
  
          return pcolor;
      }
  
      /// Core parallel plot function
      parallel_lines_handle
      axes_type::parallelplot(const std::vector<std::vector<double>> &X,
                              const std::vector<double> &colors,
                              std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          parallel_lines_handle l =
              std::make_shared<class parallel_lines>(this, X, line_spec);
          l->line_colors(colors);
          this->colormap(palette::paired());
          this->emplace_object(l);
          if (children_.size() == 1) {
              xlim({0.5, static_cast<double>(X.size()) + 0.5});
              ylim({0.0, 1.});
              this->y_axis_.visible(false);
          }
          return l;
      }
  
      /// Parallel plot - default line colors
      parallel_lines_handle
      axes_type::parallelplot(const std::vector<std::vector<double>> &X,
                              std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          return this->parallelplot(X, std::vector<double>({}), line_spec);
      }
  
      circles_handle axes_type::pie(const std::vector<double> &x,
                                    const std::vector<double> &explode,
                                    const std::vector<std::string> &labels) {
          // axes_silencer temp_silencer_{this};
  
          // pie position
          constexpr double pie_center_x = 0.0;
          constexpr double pie_center_y = 0.0;
          constexpr double pie_radius = 1.0;
          std::vector<double> x_pos(x.size(), pie_center_x);
          std::vector<double> y_pos(x.size(), pie_center_y);
          std::vector<double> radius(x.size(), pie_radius);
  
          // explode
          double real_total_x = std::accumulate(x.begin(), x.end(), 0.);
          double pie_total_x = std::max(1., real_total_x);
          if (!explode.empty()) {
              double theta_sum = 0.;
              for (size_t i = 0; i < x.size(); ++i) {
                  const double theta = theta_sum + (x[i] / pie_total_x) * pi;
                  const double explosion_radius = pie_radius * 0.1 * explode[i];
                  x_pos[i] = explosion_radius * cos(theta + pi / 2);
                  y_pos[i] = explosion_radius * sin(theta + pi / 2);
                  theta_sum += (x[i] / pie_total_x) * 2 * pi;
              }
          }
  
          // set angles
          std::vector<double> start_angle(x.size(), 0);
          std::vector<double> end_angle(x.size(), 360);
          start_angle[0] = 90;
          for (size_t i = 0; i < x.size() - 1; ++i) {
              end_angle[i] = start_angle[i] + (x[i] / pie_total_x) * 360;
              while (end_angle[i] > 360) {
                  end_angle[i] -= 360;
              }
              start_angle[i + 1] = end_angle[i];
          }
          end_angle[x.size() - 1] =
              real_total_x >= 1.0 ? 90 : (real_total_x * 360) + 90;
          while (end_angle[x.size() - 1] > 360) {
              end_angle[x.size() - 1] -= 360;
          }
  
          // pie colors
          std::vector<double> color = iota(1, x.size());
  
          // create circles in the xlim
          circles_handle pie = std::make_shared<class circles>(
              this, x_pos, y_pos, radius, start_angle, end_angle, color);
          this->emplace_object(pie);
  
          // create labels around the circles
          // convert polar coordinate to cartesian coordinates
          // t=0:0.01:2*pi;
          // x = r * cos(t);
          // y = r * sin(t);
          bool r = this->next_plot_replace();
          this->next_plot_replace(false);
          std::vector<double> label_x(x.size());
          std::vector<double> label_y(x.size());
          std::vector<std::string> pie_labels(x.size());
          double theta_sum = 0.; // [0;2pi]
          std::stringstream stream;
          for (size_t i = 0; i < x.size(); ++i) {
              if (labels.size() > i) {
                  stream << labels[i];
              } else {
                  stream << std::fixed << std::setprecision(2)
                         << x[i] * 100 / pie_total_x << "%";
              }
              pie_labels[i] = stream.str();
              stream.str("");
              double theta = theta_sum + (x[i] / pie_total_x) * pi;
              const double explosion_radius =
                  explode.size() > i ? 0.1 * explode[i] : 0;
              const double label_radius = pie_radius * (1.3 + explosion_radius);
              label_x[i] = label_radius * cos(theta + pi / 2);
              label_y[i] = label_radius * sin(theta + pi / 2);
              theta_sum += (x[i] / pie_total_x) * 2 * pi;
          }
          labels_handle ls =
              std::make_shared<class labels>(this, label_x, label_y, pie_labels);
          pie->labels(ls);
  
          // make xlim more pleasant for the pie
          if (this->x_axis().limits_mode_auto() ||
              this->y_axis().limits_mode_auto()) {
              this->axis(equal);
              this->x_axis().visible(false);
              this->y_axis().visible(false);
          }
          this->color(this->parent()->color());
          this->box(false);
  
          return pie;
      }
  
      circles_handle axes_type::pie(const std::vector<double> &x,
                                    const std::vector<std::string> &labels) {
          return this->pie(x, {}, labels);
      }
  
      line_handle axes_type::scatter(const std::vector<double> &x,
                                     const std::vector<double> &y,
                                     const std::vector<double> &sizes,
                                     const std::vector<double> &colors) {
          axes_silencer temp_silencer_{this};
          line_handle l = this->plot(x, y, "o");
          if (l->line_spec().marker_style() == line_spec::marker_style::none) {
              l->line_spec().marker_style(line_spec::marker_style::circle);
              l->line_spec().marker_face(false);
          }
          l->line_spec().line_style(line_spec::line_style::none);
          if (!sizes.empty()) {
              l->marker_size(sizes);
          }
          if (!colors.empty()) {
              l->marker_colors(colors);
          }
          return l;
      }
  
      line_handle axes_type::scatter(const std::vector<double> &x,
                                     const std::vector<double> &y, double sizes,
                                     const std::vector<double> &colors) {
          if (sizes == 0) {
              return scatter(x, y, std::vector<double>{}, colors);
          } else {
              return scatter(x, y, std::vector<double>(x.size(), sizes), colors);
          }
      }
  
      /// Core scatter3 function
      line_handle axes_type::scatter3(const std::vector<double> &x,
                                      const std::vector<double> &y,
                                      const std::vector<double> &z,
                                      const std::vector<double> &sizes,
                                      const std::vector<double> &colors,
                                      std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          line_handle l = this->plot3(x, y, z, line_spec);
          if (l->line_spec().marker_style() == line_spec::marker_style::none) {
              l->line_spec().marker_style(line_spec::marker_style::circle);
              l->line_spec().marker_face(false);
          }
          l->line_spec().line_style(line_spec::line_style::none);
          if (!sizes.empty()) {
              l->marker_size(sizes);
          }
          if (!colors.empty()) {
              l->marker_colors(colors);
          }
          return l;
      }
  
      /// Scatter3 - Default sizes and colors
      line_handle axes_type::scatter3(const std::vector<double> &x,
                                      const std::vector<double> &y,
                                      const std::vector<double> &z,
                                      std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          return this->scatter3(x, y, z, {}, {}, line_spec);
      }
  
      labels_handle
      axes_type::wordcloud(const std::vector<std::string> &words,
                           const std::vector<double> &sizes,
                           const std::vector<double> &custom_colors) {
          axes_silencer temp_silencer_{this};
  
          // set labels, their colors, and sizes give words a random position
          std::vector<std::string> labels;
          std::vector<double> colors;
          std::vector<double> float_sizes;
          double max_size = *std::max_element(sizes.begin(), sizes.end());
          constexpr double largest_word_size = 70;
          for (size_t i = 0; i < words.size(); ++i) {
              const std::string &word = words[i];
              // label is the word itself
              labels.emplace_back(word);
              // color is proportional to the frequency
              colors.emplace_back(custom_colors.empty() ? sizes[i]
                                                        : custom_colors[i]);
              // size is proportional to log(word size)
              float_sizes.emplace_back((static_cast<double>(sizes[i])) *
                                       largest_word_size / (max_size));
          }
  
          // set a position 0,0 for each label and spin until there is no overlap
          constexpr double width_factor = labels::width_factor;
          constexpr double height_factor = labels::height_factor;
          std::vector<double> x;
          std::vector<double> y;
          for (size_t i = 0; i < labels.size(); ++i) {
              x.emplace_back(0.0);
              y.emplace_back(0.0);
  
              // check if current label overlaps with previous labels
              auto overlap = [&]() {
                  const double label_width =
                      round(float_sizes[i]) * labels[i].size() * width_factor;
                  const double label_height =
                      round(float_sizes[i]) * height_factor;
                  for (size_t j = 0; j < i; ++j) {
                      const double b_label_width =
                          round(float_sizes[j]) * labels[j].size() * width_factor;
                      const double b_label_height =
                          round(float_sizes[j]) * height_factor;
                      const double axmin = x[i] - label_width / 2;
                      const double axmax = x[i] + label_width / 2;
                      const double bxmin = x[j] - b_label_width / 2;
                      const double bxmax = x[j] + b_label_width / 2;
                      const double aymin = y[i] - label_height / 2;
                      const double aymax = y[i] + label_height / 2;
                      const double bymin = y[j] - b_label_height / 2;
                      const double bymax = y[j] + b_label_height / 2;
                      if (axmin > bxmax || bxmin > axmax || aymin > bymax ||
                          bymin > aymax) {
                          // no overlap
                          continue;
                      } else {
                          // label[i] overlaps with label[j]
                          return true;
                      }
                  }
                  return false;
              };
  
              // spin the label position
              auto spin = [&]() {
                  double label_radius = sqrt(pow(x[i], 2) + pow(y[i], 2));
                  double label_theta = atan2(y[i], x[i]);
                  label_radius += 0.01;
                  label_theta += 0.1;
                  x[i] = label_radius * cos(label_theta);
                  y[i] = label_radius * sin(label_theta);
              };
  
              while (overlap()) {
                  spin();
              }
          }
  
          // attribute random size and color to each word
          // in the real word cloud
          // - position in the middle of screen
          // - rotate the word until we find a place for it
          // - create labels in the xlim
          labels_handle wordcloud = std::make_shared<class labels>(
              this, x, y, labels, colors, float_sizes);
          wordcloud->absolute_size(false);
          wordcloud->rectangles(false);
          this->color_box(true);
          this->color_box_log_scale(true);
          this->colormap(palette::default_colors_map());
          this->emplace_object(wordcloud);
  
          // make xlim more appropriate for the cloud
          this->color(this->parent()->color());
          this->box(false);
          this->x_axis().visible(false);
          this->y_axis().visible(false);
  
          return wordcloud;
      }
  
      labels_handle
      axes_type::wordcloud(const std::vector<std::string> &words,
                           const std::vector<size_t> &sizes,
                           const std::vector<double> &custom_colors) {
          return this->wordcloud(words,
                                 std::vector<double>(sizes.begin(), sizes.end()),
                                 custom_colors);
      }
  
      labels_handle
      axes_type::wordcloud(std::string_view text,
                           const std::vector<std::string> &black_list,
                           std::string_view delimiters, size_t max_cloud_size,
                           const std::vector<double> &custom_colors) {
          auto [tokens, count] =
              wordcount(text, black_list, delimiters, max_cloud_size);
          return this->wordcloud(tokens, to_vector_1d(count), custom_colors);
      }
  
      /// Core pareto function
      std::pair<bars_handle, line_handle>
      axes_type::pareto(const std::vector<double> &y,
                        const std::vector<std::string> &names, double p_threshold,
                        size_t item_threshold) {
          axes_silencer temp_silencer_{this};
          std::vector<size_t> rankings(y.size());
          rank_elements(y.begin(), y.end(), rankings.begin(),
                        std::greater<double>());
          double total_y = std::accumulate(y.begin(), y.end(), 0.);
          std::vector<std::string> x_ticklabels;
          std::vector<double> largest_y;
          std::vector<double> cummulative_y;
          for (size_t i = 0; i < rankings.size(); ++i) {
              if (names.size() > i) {
                  x_ticklabels.emplace_back(names[rankings[i]]);
              } else {
                  x_ticklabels.emplace_back(num2str(rankings[i] + 1));
              }
              largest_y.emplace_back(y[rankings[i]]);
              if (i == 0) {
                  cummulative_y.emplace_back(largest_y[i]);
              } else {
                  cummulative_y.emplace_back(cummulative_y[i - 1] + largest_y[i]);
              }
              if (cummulative_y[i] >= p_threshold * total_y ||
                  i + 1 >= item_threshold) {
                  if (i < rankings.size() - 1) {
                      cummulative_y.emplace_back(cummulative_y[i] +
                                                 y[rankings[i + 1]]);
                  }
                  break;
              }
          }
          cummulative_y = transform(
              cummulative_y, [total_y](double x) { return x * 100. / total_y; });
  
          // plot line
          line_handle l = plot(cummulative_y);
          l->line_width(1.);
          this->y2_axis_.visible(true);
          l->use_y2(true);
  
          // plot bars
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
          bars_handle b = this->bar(largest_y);
          b->face_color(l->color());
          this->next_plot_replace(p2);
  
          // set range
          this->y_axis().limits({0, total_y});
  
          // set ticks
          this->x_axis().tick_values(iota(1, largest_y.size()));
          this->x_axis().ticklabels(x_ticklabels);
          this->x_axis().limits({0.5, double(x_ticklabels.size()) + 0.5});
  
          // set y2 axis
          this->y2_axis().visible(true);
          this->y2_axis().limits({0, inf});
          this->y2_axis().tick_label_format("%g%%");
  
          this->box(true);
  
          return std::make_pair(b, l);
      }
  
      /// Core plot function
      line_handle axes_type::stem(const std::vector<double> &x,
                                  const std::vector<double> &y,
                                  std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          line_handle l = this->plot(x, y, line_spec);
          l->impulse(true);
          this->x_axis().zero_axis(true);
          return l;
      }
  
      /// Basic plot function with automatic x
      line_handle axes_type::stem(const std::vector<double> &y,
                                  std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          line_handle l = this->stem(iota(1, y.size()), y, line_spec);
          l->impulse(true);
          this->x_axis().zero_axis(true);
  
          return l;
      }
  
      /// Stem - Plot lists of lists
      std::vector<line_handle>
      axes_type::stem(const std::vector<double> &x,
                      const std::vector<std::vector<double>> &Y,
                      std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          std::vector<line_handle> hs = this->plot(x, Y, line_spec);
          for (auto &h : hs) {
              h->impulse(true);
          }
          this->x_axis().zero_axis(true);
  
          return hs;
      }
  
      /// Plot lists of lists with automatic x
      std::vector<line_handle>
      axes_type::stem(const std::vector<std::vector<double>> &Y,
                      std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          std::vector<line_handle> hs = this->plot(Y, line_spec);
          for (auto &h : hs) {
              h->impulse(true);
          }
          this->x_axis().zero_axis(true);
  
          return hs;
      }
  
      /// Stem 3d - Core function
      line_handle axes_type::stem3(const std::vector<double> &x,
                                   const std::vector<double> &y,
                                   const std::vector<double> &z,
                                   std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          line_handle l = this->plot3(x, y, z, line_spec);
          l->impulse(true);
          return l;
      }
  
      /// 3d-stem lists of Xs and Ys
      std::vector<line_handle>
      axes_type::stem3(const std::vector<std::vector<double>> &X,
                       const std::vector<std::vector<double>> &Y,
                       const std::vector<double> &z,
                       std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          auto hs = this->plot3(X, Y, z, line_spec);
          for (auto &h : hs) {
              h->impulse(true);
          }
          return hs;
      }
  
      /// Stem 3d - Automatic x and y
      line_handle axes_type::stem3(const std::vector<double> &z,
                                   std::string_view line_spec) {
          return this->stem3(iota(1, z.size()), std::vector<double>(z.size(), 1.),
                             z, line_spec);
      }
  
      /// Stem 3d - Automatic x and y - Many Zs
      line_handle axes_type::stem3(const std::vector<std::vector<double>> &Z,
                                   std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          std::vector<double> x;
          std::vector<double> y;
          std::vector<double> z;
          for (size_t i = 0; i < Z.size(); ++i) {
              auto tmp_x = iota(1, Z[i].size());
              x.insert(x.end(), tmp_x.begin(), tmp_x.end());
              std::vector<double> tmp_y(Z[i].size(), i);
              y.insert(y.end(), tmp_y.begin(), tmp_y.end());
              z.insert(z.end(), Z[i].begin(), Z[i].end());
          }
          return this->stem3(x, y, z, line_spec);
      }
  
      /// Core geoplot function
      circles_handle axes_type::geobubble(const std::vector<double> &latitude,
                                          const std::vector<double> &longitude,
                                          const std::vector<double> &sizes,
                                          const std::vector<double> &colors) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
  
          // create map if we need
          this->geoplot();
  
          // normalize sizes
          constexpr double max_bubble_size = 10.;
          std::vector<double> normalized_sizes(sizes.begin(), sizes.end());
          auto mit = std::max_element(sizes.begin(), sizes.end());
          if (mit != sizes.end()) {
              double max_size = *mit;
              for (auto &size : normalized_sizes) {
                  size /= max_size;
                  size *= max_bubble_size;
              }
          }
  
          this->next_plot_replace(false);
          circles_handle h = std::make_shared<circles>(
              this, longitude, latitude, normalized_sizes, std::vector<double>(),
              std::vector<double>(), colors);
          h->line_width(1.);
          auto fc = this->get_color_and_bump();
          const double alpha = fc[0];
          const double opacity = (1 - fc[0]);
          const double new_opacity = opacity * 0.7;
          const double new_alpha = 1 - new_opacity;
          fc[0] = static_cast<float>(new_alpha);
          h->face_color(fc);
          h->line_color(to_array("white"));
          this->emplace_object(h);
  
          this->next_plot_replace(p2);
          return h;
      }
  
      /// Core geoplot function
      line_handle axes_type::geodensityplot(const std::vector<double> &latitude,
                                            const std::vector<double> &longitude,
                                            const std::vector<double> &weights) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
  
          // create map if we need
          this->geoplot();
  
          // weighted bin count
          auto [min_x, max_x] =
              std::minmax_element(longitude.begin(), longitude.end());
          auto x_edges = histogram::bin_picker(*min_x, *max_x, 200, 0);
  
          auto [min_y, max_y] =
              std::minmax_element(latitude.begin(), latitude.end());
          auto y_edges = histogram::bin_picker(*min_y, *max_y, 200, 0);
  
          std::vector<std::vector<double>> bin_counts(
              x_edges.size() - 1, std::vector<double>(y_edges.size() - 1, 0.));
          for (size_t i = 0; i < longitude.size(); ++i) {
              double xi = longitude[i];
              double yi = latitude[i];
              double wi = weights.size() > i ? weights[i] : 1.;
              // find first edge that does not compare less than v
              auto itx = std::lower_bound(x_edges.begin(), x_edges.end(), xi);
              bool xout_of_range = itx == x_edges.begin() || itx == x_edges.end();
              if (!xout_of_range) {
                  auto ity = std::lower_bound(y_edges.begin(), y_edges.end(), yi);
                  bool yout_of_range =
                      ity == y_edges.begin() || ity == y_edges.end();
                  if (!yout_of_range) {
                      bin_counts[itx - x_edges.begin() - 1]
                                [ity - y_edges.begin() - 1] += wi;
                  }
              }
          }
  
          // create bin positions and their sizes
          std::vector<double> bin_x;
          std::vector<double> bin_y;
          std::vector<double> colors;
          for (size_t i = 0; i < x_edges.size() - 1; ++i) {
              for (size_t j = 0; j < y_edges.size() - 1; ++j) {
                  if (bin_counts[i][j] != 0.0) {
                      bin_x.emplace_back((x_edges[i] + x_edges[i + 1]) / 2.);
                      bin_y.emplace_back((y_edges[j] + y_edges[j + 1]) / 2.);
                      colors.emplace_back(bin_counts[i][j]);
                  }
              }
          }
  
          this->next_plot_replace(false);
          line_handle s = this->scatter(bin_x, bin_y, {}, colors);
          s->marker_face(true);
          s->marker_size(1.5);
          this->color_box(true);
  
          this->next_plot_replace(p2);
          return s;
      }
  
      /// Plot world map on the axes
      line_handle axes_type::geoplot() {
          axes_silencer temp_silencer_{this};
  
          for (const auto &child : this->children()) {
              if (child->tag() == "map") {
                  return std::dynamic_pointer_cast<class line>(child);
              }
          }
  
          auto [x, y] = world_map_110m();
          line_handle a = this->plot(x, y);
          a->tag("map");
          color_array land_color = {0, 0.9294, 0.9294, 0.9294};
          a->color(land_color);
          a->fill(true);
          color_array bg = {0, 0.7882, 0.7882, 0.7882};
          this->color(bg);
  
          this->x_axis().geographic(true);
          this->x_axis().tick_label_format("%Dº%E");
          this->x_axis().label("Longitude");
          this->y_axis().geographic(true);
          this->y_axis().tick_label_format("%Dº%N");
          this->y_axis().label("Latitude");
  
          return a;
      }
  
      /// Core geoplot function - Plot lines on world map
      line_handle axes_type::geoplot(const std::vector<double> &latitude,
                                     const std::vector<double> &longitude,
                                     std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
          this->geoplot();
          this->next_plot_replace(false);
          line_handle l = this->plot(longitude, latitude, line_spec);
          l->line_width(1.0);
          this->next_plot_replace(p2);
          return l;
      }
  
      /// Adjust limits and reload that region of the world map in the proper
      /// resolution
      void axes_type::geolimits(const std::array<double, 2> &latitude,
                                const std::array<double, 2> &longitude) {
          axes_silencer temp_silencer_{this};
  
          // look for map object in the xlim
          line_handle map = nullptr;
          for (const axes_object_handle &c : this->children()) {
              if (c->tag() == "map") {
                  map = std::dynamic_pointer_cast<class line>(c);
              }
          }
  
          // If there is no map in the axes
          if (!map) {
              map = this->geoplot();
          }
  
          // if we found the map
          if (map) {
              // calculate the perimeter of the new limits
              double latitude_kms = 111.12 * std::abs(latitude[1] - latitude[0]);
              double longitude_kms =
                  111.12 * std::abs(longitude[1] - longitude[0]);
              double w_pixels = this->width() * this->parent()->width();
              double h_pixels = this->height() * this->parent()->height();
              double w_km_per_pixel = longitude_kms / w_pixels;
              double h_km_per_pixel = latitude_kms / h_pixels;
              double min_km_per_pixel = std::min(w_km_per_pixel, h_km_per_pixel);
  
              auto &[map_x, map_y] = (min_km_per_pixel <= 10) ? world_map_10m()
                                     : (min_km_per_pixel <= 50)
                                         ? world_map_50m()
                                         : world_map_110m();
              std::vector<double> limits_map_x;
              std::vector<double> limits_map_y;
  
              // for each submap
              for (size_t i = 0; i < map_x.size(); ++i) {
                  // copy next submap
                  std::vector<double> submap_x;
                  std::vector<double> submap_y;
                  while (i < map_x.size() && std::isfinite(map_x[i]) &&
                         std::isfinite(map_y[i])) {
                      submap_x.emplace_back(map_x[i]);
                      submap_y.emplace_back(map_y[i]);
                      ++i;
                  }
  
                  // fn to check if a point is inside the limits we want
                  auto inside_the_map = [&](double x, double y) {
                      return std::isfinite(x) && std::isfinite(y) &&
                             x <= longitude[1] && x >= longitude[0] &&
                             y <= latitude[1] && y >= latitude[0];
                  };
  
                  // fn to check in what quadrant a point falls around our limits
                  auto quadrant = [&](double x, double y) {
                      if (y > latitude[1]) {
                          if (x < longitude[0]) {
                              return 1;
                          } else if (x > longitude[1]) {
                              return 3;
                          } else {
                              return 2;
                          }
                      } else if (y < latitude[0]) {
                          if (x < longitude[0]) {
                              return 7;
                          } else if (x > longitude[1]) {
                              return 9;
                          } else {
                              return 8;
                          }
                      } else {
                          if (x < longitude[0]) {
                              return 4;
                          } else if (x > longitude[1]) {
                              return 6;
                          } else {
                              return 5;
                          }
                      }
                  };
  
                  // check if submap has points inside the limits
                  bool add_submap = false;
                  for (size_t j = 0; j < submap_x.size(); ++j) {
                      if (inside_the_map(submap_x[j], submap_y[j])) {
                          add_submap = true;
                          break;
                      }
                  }
  
                  // copy submap if it has points we need
                  if (add_submap) {
                      int last_q = quadrant(submap_x[0], submap_y[0]);
                      int q = last_q;
                      // for each point
                      for (size_t j = 0; j < submap_x.size(); ++j) {
                          // if point is inside the limit
                          if (inside_the_map(submap_x[j], submap_y[j])) {
                              // copy point
                              q = 5;
                              limits_map_x.emplace_back(submap_x[j]);
                              limits_map_y.emplace_back(submap_y[j]);
                          } else {
                              // else, if changed the quadrant around the map
                              q = quadrant(submap_x[j], submap_y[j]);
                              if (q != last_q) {
                                  // copy point
                                  limits_map_x.emplace_back(submap_x[j]);
                                  limits_map_y.emplace_back(submap_y[j]);
                              }
                          }
                          last_q = q;
                      }
                      // new submap divider
                      if (i != map_x.size() - 1) {
                          limits_map_x.emplace_back(NaN);
                          limits_map_y.emplace_back(NaN);
                      }
                  }
              }
  
              map->x_data(limits_map_x);
              map->y_data(limits_map_y);
          }
  
          this->x_axis().limits(longitude);
          this->y_axis().limits(latitude);
      }
  
      void axes_type::geolimits(double latitude_x, double latitude_y,
                                double longitude_x, double longitude_y) {
          geolimits(std::array<double, 2>{latitude_x, latitude_y},
                    std::array<double, 2>{longitude_x, longitude_y});
      }
  
      void axes_type::geolimits(const std::vector<double> &latitude,
                                const std::vector<double> &longitude) {
          this->geolimits(to_array<2>(latitude), to_array<2>(longitude));
      }
  
      /// Core geoscatter function
      line_handle axes_type::geoscatter(const std::vector<double> &latitude,
                                        const std::vector<double> &longitude,
                                        const std::vector<double> &sizes,
                                        const std::vector<double> &colors) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
  
          // create map if we need
          this->geoplot();
  
          this->next_plot_replace(false);
          line_handle l = this->scatter(longitude, latitude, sizes, colors);
  
          this->next_plot_replace(p2);
          return l;
      }
  
      void axes_type::axis(const std::array<double, 4> &limits_x_y) {
          this->x_axis().limits({limits_x_y[0], limits_x_y[1]});
          this->x_axis().limits_mode_auto(false);
          this->y_axis().limits({limits_x_y[2], limits_x_y[3]});
          this->y_axis().limits_mode_auto(false);
          this->touch();
      }
  
      void axes_type::axis(bool v) {
          this->visible(v);
          if (!v) {
              this->title_visible(true);
          }
      }
  
      /// \brief Make x/y/z-axis limits automatic
      void axes_type::axis(keyword_automatic_type automatic) {
          this->limits_mode(automatic);
      }
  
      /// \brief Make x/y/z-axis limits manual
      void axes_type::axis(keyword_manual_type manual) {
          this->limits_mode(manual);
      }
  
      /// \brief Reverse y-axis
      void axes_type::axis(keyword_ij_type) { this->y_axis().reverse(true); }
  
      /// \brief Make x/y-axis ratio = 1 by changing the axis limits
      void axes_type::axis(keyword_equal_type) {
          if (this->children().empty()) {
              return;
          }
          float w = this->width() * this->parent()->width();
          float h = this->height() * this->parent()->height();
          auto [xmin, xmax, ymin, ymax] = this->child_limits();
          double xrange = xmax - xmin;
          double yrange = ymax - ymin;
          // new ranges that could make it proportional
          // new xrange = w / factor that makes current xrange fit in there
          // w / factor == xrange -> factor = w / xrange
          double new_xrange = w * (yrange / h);
          double new_yrange = h * (xrange / w);
          if (new_xrange >= xrange) {
              double size_increase = (new_xrange - xrange) / 2;
              this->x_axis().limits({xmin - size_increase, xmax + size_increase});
              this->y_axis().limits({ymin, ymax});
          } else if (new_yrange >= yrange) {
              double size_increase = (new_yrange - yrange) / 2;
              this->x_axis().limits({xmin, xmax});
              this->y_axis().limits({ymin - size_increase, ymax + size_increase});
          }
      }
  
      /// \brief Make x/y/z-axis limits tight from x/y/z-min to x/y/z-max
      void axes_type::axis(keyword_tight_type) {
          if (this->children().empty()) {
              return;
          }
          auto [xmin, xmax, ymin, ymax] = this->child_limits();
          this->x_axis().limits({xmin, xmax});
          this->y_axis().limits({ymin, ymax});
          if (this->is_3d() && !this->is_3d_map()) {
              auto zminit = std::min_element(
                  this->children().begin(), this->children().end(),
                  [](axes_object_handle a, axes_object_handle b) {
                      return a->zmin() < b->zmin();
                  });
              auto zmaxit = std::max_element(
                  this->children().begin(), this->children().end(),
                  [](axes_object_handle a, axes_object_handle b) {
                      return a->zmax() < b->zmax();
                  });
              this->z_axis().limits({(*zminit)->zmin(), (*zmaxit)->zmax()});
          }
      }
  
      /// \brief Make x/y-axis square with ratio 1 by changing the axes size
      void axes_type::axis(keyword_square_type) {
          axis(equal);
          double w = this->width() * this->parent()->width();
          double h = this->height() * this->parent()->height();
          // make w = h
          if (w > h) {
              double old_width = this->width();
              double new_width = h / this->parent()->width();
              this->width(static_cast<float>(new_width));
              this->x_origin(this->x_origin() +
                             static_cast<float>(old_width - new_width) / 2.f);
          } else if (h > w) {
              double old_height = this->height();
              double new_height = h / this->parent()->height();
              this->height(static_cast<float>(new_height));
              this->y_origin(this->y_origin() + static_cast<float>(old_height - new_height) / 2.f);
          }
      }
  
      /// Compass plot function
      vectors_handle axes_type::compass(const std::vector<double> &x,
                                        const std::vector<double> &y,
                                        std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          auto theta = transform(
              x, y, [](double x, double y) { return vector_radians(x, y); });
          auto rho = transform(
              x, y, [](double x, double y) { return vector_magnitude(x, y); });
  
          vectors_handle l =
              std::make_shared<class vectors>(this, theta, rho, line_spec);
          this->emplace_object(l);
  
          l->polar(true);
  
          this->axis(equal);
          this->x_axis().visible(false);
          this->y_axis().visible(false);
          this->r_axis().visible(false);
          this->r_axis().tick_length(0.);
          this->t_axis().visible(true);
          this->t_axis().tick_length(0);
  
          return l;
      }
  
      /// Ezpolar - Plot function on polar plot
      string_function_handle axes_type::ezpolar(std::string_view equation,
                                                std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          string_function_handle obj =
              std::make_shared<class string_function>(this, equation, line_spec);
  
          this->emplace_object(obj);
          obj->polar(true);
  
          this->axis(equal);
          this->x_axis().visible(false);
          this->y_axis().visible(false);
          this->r_axis().visible(false);
          this->r_axis().tick_length(0.);
          this->t_axis().visible(true);
          this->t_axis().tick_length(0);
  
          return obj;
      }
  
      std::vector<string_function_handle>
      axes_type::ezpolar(std::vector<std::string> equations,
                         std::vector<std::string> line_specs) {
          std::vector<string_function_handle> res;
          auto it_line_specs = line_specs.begin();
          for (const auto &equation : equations) {
              if (it_line_specs != line_specs.end()) {
                  res.emplace_back(this->ezpolar(equation, *it_line_specs));
              } else {
                  res.emplace_back(this->ezpolar(equation));
              }
              ++it_line_specs;
          }
          return res;
      }
  
      function_line_handle
      axes_type::ezpolar(function_line::function_type equation,
                         const std::array<double, 2> &t_range,
                         std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          function_line_handle obj = std::make_shared<class function_line>(
              this, equation, t_range, line_spec);
          this->emplace_object(obj);
          obj->polar(true);
  
          this->axis(equal);
          this->x_axis().visible(false);
          this->y_axis().visible(false);
          this->r_axis().visible(false);
          this->r_axis().tick_length(0.);
          this->t_axis().visible(true);
          this->t_axis().tick_length(0);
  
          return obj;
      }
  
      function_line_handle
      axes_type::ezpolar(function_line::function_type equation,
                         std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          function_line_handle obj = std::make_shared<class function_line>(
              this, equation, std::array<double, 2>({0, 2 * pi}), line_spec);
          this->emplace_object(obj);
          obj->polar(true);
          this->axis(equal);
          this->x_axis().visible(false);
          this->y_axis().visible(false);
          this->r_axis().visible(false);
          this->r_axis().tick_length(0.);
          this->t_axis().visible(true);
          this->t_axis().tick_length(0);
  
          return obj;
      }
  
      function_line_handle
      axes_type::ezpolar(function_line::function_type equation,
                         std::vector<double> x_range,
                         std::string_view line_spec) {
          return this->ezpolar(equation, to_array<2>(x_range), line_spec);
      }
  
      function_line_handle
      axes_type::ezpolar(function_line::function_type function_x,
                         function_line::function_type function_y,
                         const std::array<double, 2> &t_range,
                         std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          function_line_handle obj = std::make_shared<class function_line>(
              this, function_x, function_y, t_range, line_spec);
          this->emplace_object(obj);
          obj->polar(true);
          this->axis(equal);
          this->x_axis().visible(false);
          this->y_axis().visible(false);
          this->r_axis().visible(false);
          this->r_axis().tick_length(0.);
          this->t_axis().visible(true);
          this->t_axis().tick_length(0);
  
          return obj;
      }
  
      function_line_handle
      axes_type::ezpolar(function_line::function_type function_x,
                         function_line::function_type function_y,
                         std::vector<double> t_range,
                         std::string_view line_spec) {
          return this->ezpolar(function_x, function_y, to_array<2>(t_range),
                               line_spec);
      }
  
      std::vector<function_line_handle>
      axes_type::ezpolar(std::vector<function_line::function_type> equations,
                         std::array<double, 2> x_range,
                         std::vector<std::string> line_specs) {
          axes_silencer temp_silencer_{this};
          std::vector<function_line_handle> res;
          auto it_line_specs = line_specs.begin();
          for (const auto &equation : equations) {
              if (it_line_specs != line_specs.end()) {
                  res.emplace_back(
                      this->ezpolar(equation, x_range, *it_line_specs));
              } else {
                  res.emplace_back(this->ezpolar(equation, x_range));
              }
              ++it_line_specs;
          }
          return res;
      }
  
      std::vector<function_line_handle>
      axes_type::ezpolar(std::vector<function_line::function_type> equations,
                         std::vector<double> x_range,
                         std::vector<std::string> line_specs) {
          return this->ezpolar(equations, to_array<2>(x_range), line_specs);
      }
  
      /// Polar histogram - Core function
      histogram_handle axes_type::polarhistogram(const std::vector<double> &theta,
                                                 size_t nbins) {
          axes_silencer temp_silencer_{this};
  
          auto reflect_theta =
              transform(theta, [](double t) { return to_positive_radian(t); });
          auto edges = linspace(0., 2 * pi, nbins + 1);
          histogram_handle h = this->hist(reflect_theta, edges);
          h->polar(true);
          this->axis(equal);
  
          this->x_axis().visible(false);
          this->y_axis().visible(false);
          this->r_axis().visible(false);
          this->r_axis().tick_length(0.);
          this->t_axis().visible(true);
          this->t_axis().tick_length(0);
  
          return h;
      }
  
      /// Polar plot - Core function
      line_handle axes_type::polarplot(const std::vector<double> &theta,
                                       const std::vector<double> &rho,
                                       std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          line_handle l = this->plot(theta, rho, line_spec);
          l->polar(true);
          this->axis(equal);
  
          this->x_axis().visible(false);
          this->y_axis().visible(false);
          this->r_axis().visible(false);
          this->r_axis().tick_length(0.);
          this->t_axis().visible(true);
          this->t_axis().tick_length(0);
  
          return l;
      }
  
      /// Polar plot - Automatic theta
      line_handle axes_type::polarplot(const std::vector<double> &rho,
                                       std::string_view line_spec) {
          return this->polarplot(linspace(0, 2 * pi, rho.size()), rho, line_spec);
      }
  
      /// Polar plot - Complex numbers
      line_handle axes_type::polarplot(const std::vector<std::complex<double>> &z,
                                       std::string_view line_spec) {
          std::vector<double> theta;
          std::vector<double> rho;
          for (size_t i = 0; i < z.size(); ++i) {
              rho.emplace_back(sqrt(pow(z[i].real(), 2) + pow(z[i].imag(), 2)));
              theta.emplace_back(atan2(z[i].imag(), z[i].real()));
          }
          return this->polarplot(theta, rho, line_spec);
      }
  
      /// Polar scatter - Core function
      line_handle axes_type::polarscatter(const std::vector<double> &theta,
                                          const std::vector<double> &rho,
                                          const std::vector<double> &sizes,
                                          const std::vector<double> &colors,
                                          std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          line_handle l = this->scatter(theta, rho, sizes, colors);
          l->line_style(line_spec);
          l->polar(true);
          this->axis(equal);
  
          this->x_axis().visible(false);
          this->y_axis().visible(false);
          this->r_axis().visible(false);
          this->r_axis().tick_length(0.);
          this->t_axis().visible(true);
          this->t_axis().tick_length(0);
  
          return l;
      }
  
      /// Polar scatter - Single size - Default colors
      line_handle axes_type::polarscatter(const std::vector<double> &theta,
                                          const std::vector<double> &rho,
                                          double size,
                                          std::string_view line_spec) {
          return this->polarscatter(theta, rho, std::vector(theta.size(), size),
                                    {}, line_spec);
      }
  
      /// Polar scatter - Default size and colors
      line_handle axes_type::polarscatter(const std::vector<double> &theta,
                                          const std::vector<double> &rho,
                                          std::string_view line_spec) {
          return this->polarscatter(theta, rho, {}, {}, line_spec);
      }
  
      /// Contour - Core function - Manual levels
      contours_handle
      axes_type::contour(const std::vector<std::vector<double>> &X,
                         const std::vector<std::vector<double>> &Y,
                         const std::vector<std::vector<double>> &Z,
                         std::vector<double> levels, std::string_view line_spec,
                         size_t n_levels) {
          axes_silencer temp_silencer_{this};
  
          contours_handle l =
              std::make_shared<class contours>(this, X, Y, Z, line_spec);
          if (n_levels) {
              l->n_levels(n_levels);
          }
          if (!levels.empty()) {
              l->levels(levels);
          }
          this->x_axis().limits({l->xmin(), l->xmax()});
          this->y_axis().limits({l->ymin(), l->ymax()});
          this->emplace_object(l);
          this->color_box(true);
  
          return l;
      }
  
      /// Contour - Manual number of levels (or 0 for automatic number of levels)
      contours_handle
      axes_type::contour(const std::vector<std::vector<double>> &X,
                         const std::vector<std::vector<double>> &Y,
                         const std::vector<std::vector<double>> &Z,
                         size_t n_levels, std::string_view line_spec) {
          return this->contour(X, Y, Z, {}, line_spec, n_levels);
      }
  
      /// Contour - Automatic levels and number of levels
      contours_handle
      axes_type::contour(const std::vector<std::vector<double>> &X,
                         const std::vector<std::vector<double>> &Y,
                         const std::vector<std::vector<double>> &Z,
                         std::string_view line_spec) {
          return this->contour(X, Y, Z, 0, line_spec);
      }
  
      /// Contour filled - Manual levels
      contours_handle
      axes_type::contourf(const std::vector<std::vector<double>> &X,
                          const std::vector<std::vector<double>> &Y,
                          const std::vector<std::vector<double>> &Z,
                          std::vector<double> levels,
                          std::string_view line_spec, size_t n_levels) {
          axes_silencer temp_silencer_{this};
  
          contours_handle l = this->contour(X, Y, Z, levels, line_spec, n_levels);
          l->filled(true);
          l->line_style("k");
          this->emplace_object(l);
  
          return l;
      }
  
      /// Contour filled - Manual number of levels
      contours_handle
      axes_type::contourf(const std::vector<std::vector<double>> &X,
                          const std::vector<std::vector<double>> &Y,
                          const std::vector<std::vector<double>> &Z,
                          size_t n_levels, std::string_view line_spec) {
          return this->contourf(X, Y, Z, {}, line_spec, n_levels);
      }
  
      /// Contour filled - Automatic number of levels
      contours_handle
      axes_type::contourf(const std::vector<std::vector<double>> &X,
                          const std::vector<std::vector<double>> &Y,
                          const std::vector<std::vector<double>> &Z,
                          std::string_view line_spec) {
          return this->contourf(X, Y, Z, 0, line_spec);
      }
  
      using fcontour_function_type = std::function<double(double, double)>;
  
      /// Lambda function contour - Manual levels (or empty list for automatic) /
      /// Manual number of levels (or 0 for automatic)
      contours_handle axes_type::fcontour(fcontour_function_type fn,
                                          const std::array<double, 4> &xy_range,
                                          std::vector<double> levels,
                                          std::string_view line_spec,
                                          size_t n_levels) {
          axes_silencer temp_silencer_{this};
  
          std::vector<double> x = linspace(xy_range[0], xy_range[1]);
          std::vector<double> y = linspace(xy_range[2], xy_range[3]);
          auto [X, Y] = meshgrid(x, y);
          vector_2d Z(y.size(), vector_1d(x.size(), 0.));
          for (size_t i = 0; i < y.size(); ++i) {
              for (size_t j = 0; j < x.size(); ++j) {
                  Z[i][j] = fn(X[i][j], Y[i][j]);
              }
          }
  
          contours_handle l = this->contour(X, Y, Z, levels, line_spec, n_levels);
  
          return l;
      }
  
      /// Lambda function contour - Manual number of levels (default = 9)
      contours_handle axes_type::fcontour(fcontour_function_type fn,
                                          size_t n_levels,
                                          std::string_view line_spec) {
          return this->fcontour(fn, std::array<double, 4>{-5, +5, -5, +5},
                                std::vector<double>{}, line_spec, n_levels);
      }
  
      /// Lambda function contour - Automatic number of levels and levels
      contours_handle axes_type::fcontour(fcontour_function_type fn,
                                          std::string_view line_spec) {
          return this->fcontour(fn, 9, line_spec);
      }
  
      /// Feather - Core function
      vectors_handle axes_type::feather(const std::vector<double> &u,
                                        const std::vector<double> &v,
                                        std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          vectors_handle l =
              std::make_shared<class vectors>(this, u, v, line_spec);
          this->emplace_object(l);
  
          if (!l->line_spec().user_color()) {
              l->line_spec().color(this->get_color_and_bump());
          }
  
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
          std::vector<double> base_line(u.size(), 0);
          auto b = this->plot(base_line, line_spec);
          b->line_spec().color(l->line_spec().color());
          this->next_plot_replace(p2);
  
          return l;
      }
  
      /// Quiver - Core function
      vectors_handle axes_type::quiver(const std::vector<double> &x,
                                       const std::vector<double> &y,
                                       const std::vector<double> &u,
                                       const std::vector<double> &v, double scale,
                                       std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          auto x_copy = x;
          std::sort(x_copy.begin(), x_copy.end());
          x_copy.resize(std::distance(x_copy.begin(),
                                      std::unique(x_copy.begin(), x_copy.end())));
          std::vector<double> x_diff(x_copy.size());
          std::adjacent_difference(x_copy.begin(), x_copy.end(), x_diff.begin());
          auto xdiffmin = std::min_element(x_diff.begin() + 1, x_diff.end());
          double u_max = *std::max_element(u.begin(), u.end());
          double x_max = xdiffmin != x_diff.end() ? *xdiffmin : u_max;
  
          auto y_copy = y;
          std::sort(y_copy.begin(), y_copy.end());
          y_copy.resize(std::distance(y_copy.begin(),
                                      std::unique(y_copy.begin(), y_copy.end())));
          std::vector<double> y_diff(y_copy.size());
          std::adjacent_difference(y_copy.begin(), y_copy.end(), y_diff.begin());
          auto ydiffmin = std::min_element(y_diff.begin() + 1, y_diff.end());
          double v_max = *std::max_element(v.begin(), v.end());
          double y_max = ydiffmin != y_diff.end() ? *ydiffmin : v_max;
  
          auto u_scaled =
              (scale != 0.)
                  ? transform(
                        u, [&](double u) { return (u / u_max) * scale * x_max; })
                  : u;
          auto v_scaled =
              (scale != 0.)
                  ? transform(
                        v, [&](double v) { return (v / v_max) * scale * y_max; })
                  : v;
  
          vectors_handle l = std::make_shared<class vectors>(this, x, y, u_scaled,
                                                             v_scaled, line_spec);
          this->emplace_object(l);
  
          return l;
      }
  
      /// Quiver - 2d x,y,u,v
      vectors_handle axes_type::quiver(const std::vector<std::vector<double>> &x,
                                       const std::vector<std::vector<double>> &y,
                                       const std::vector<std::vector<double>> &u,
                                       const std::vector<std::vector<double>> &v,
                                       double scale,
                                       std::string_view line_spec) {
          return this->quiver(flatten(x), flatten(y), flatten(u), flatten(v),
                              scale, line_spec);
      }
  
      /// Quiver 3d - Core function
      vectors_handle axes_type::quiver3(
          const std::vector<double> &x, const std::vector<double> &y,
          const std::vector<double> &z, const std::vector<double> &u,
          const std::vector<double> &v, const std::vector<double> &w,
          double scale, std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          auto x_copy = x;
          std::sort(x_copy.begin(), x_copy.end());
          x_copy.resize(std::distance(x_copy.begin(),
                                      std::unique(x_copy.begin(), x_copy.end())));
          std::vector<double> x_diff(x_copy.size());
          std::adjacent_difference(x_copy.begin(), x_copy.end(), x_diff.begin());
          auto xdiffmin = std::min_element(x_diff.begin() + 1, x_diff.end());
          double u_max = *std::max_element(u.begin(), u.end());
          double x_max = xdiffmin != x_diff.end() ? *xdiffmin : u_max;
  
          auto y_copy = y;
          std::sort(y_copy.begin(), y_copy.end());
          y_copy.resize(std::distance(y_copy.begin(),
                                      std::unique(y_copy.begin(), y_copy.end())));
          std::vector<double> y_diff(y_copy.size());
          std::adjacent_difference(y_copy.begin(), y_copy.end(), y_diff.begin());
          auto ydiffmin = std::min_element(y_diff.begin() + 1, y_diff.end());
          double v_max = *std::max_element(v.begin(), v.end());
          double y_max = ydiffmin != y_diff.end() ? *ydiffmin : v_max;
  
          auto z_copy = y;
          std::sort(z_copy.begin(), z_copy.end());
          z_copy.resize(std::distance(z_copy.begin(),
                                      std::unique(z_copy.begin(), z_copy.end())));
          std::vector<double> z_diff(z_copy.size());
          std::adjacent_difference(z_copy.begin(), z_copy.end(), z_diff.begin());
          auto zdiffmin = std::min_element(z_diff.begin() + 1, z_diff.end());
          double w_max = *std::max_element(v.begin(), v.end());
          double z_max = ydiffmin != z_diff.end() ? *ydiffmin : w_max;
  
          auto u_scaled =
              (scale != 0.)
                  ? transform(
                        u, [&](double u) { return (u / u_max) * scale * x_max; })
                  : u;
          auto v_scaled =
              (scale != 0.)
                  ? transform(
                        v, [&](double v) { return (v / v_max) * scale * y_max; })
                  : v;
          auto w_scaled =
              (scale != 0.)
                  ? transform(
                        w, [&](double w) { return (w / w_max) * scale * z_max; })
                  : w;
  
          vectors_handle l = std::make_shared<class vectors>(
              this, x, y, z, u_scaled, v_scaled, w_scaled, line_spec);
          this->emplace_object(l);
  
          return l;
      }
  
      /// Quiver 3d - 2d vectors
      vectors_handle axes_type::quiver3(const std::vector<std::vector<double>> &x,
                                        const std::vector<std::vector<double>> &y,
                                        const std::vector<std::vector<double>> &z,
                                        const std::vector<std::vector<double>> &u,
                                        const std::vector<std::vector<double>> &v,
                                        const std::vector<std::vector<double>> &w,
                                        double scale,
                                        std::string_view line_spec) {
          return this->quiver3(flatten(x), flatten(y), flatten(z), flatten(u),
                               flatten(v), flatten(w), scale, line_spec);
      }
  
      /// Quiver 3d - Automatic x and y - 2d vectors
      vectors_handle axes_type::quiver3(const std::vector<std::vector<double>> &z,
                                        const std::vector<std::vector<double>> &u,
                                        const std::vector<std::vector<double>> &v,
                                        const std::vector<std::vector<double>> &w,
                                        double scale,
                                        std::string_view line_spec) {
          auto [n, m] = size(z);
          vector_1d x = iota(1, m);
          vector_1d y = iota(1, n);
          auto [xx, yy] = meshgrid(x, y);
          return this->quiver3(xx, yy, z, u, v, w, scale, line_spec);
      }
  
      /// Fence - Core function
      surface_handle axes_type::fence(const std::vector<std::vector<double>> &X,
                                      const std::vector<std::vector<double>> &Y,
                                      const std::vector<std::vector<double>> &Z,
                                      const std::vector<double> &c) {
          axes_silencer temp_silencer_{this};
  
          surface_handle l =
              this->mesh(X, Y, Z, !c.empty() ? vector_2d{c} : vector_2d{});
          l->fences(true);
  
          return l;
      }
  
      // z = f(x,y)
      using fmesh_function_type = std::function<double(double, double)>;
  
      /// Lambda function mesh - Core function
      surface_handle axes_type::fmesh(fcontour_function_type fn,
                                      const std::array<double, 2> &x_range,
                                      const std::array<double, 2> &y_range,
                                      double mesh_density) {
          axes_silencer temp_silencer_{this};
  
          std::vector<double> x =
              linspace(x_range[0], x_range[1], static_cast<size_t>(mesh_density));
          std::vector<double> y =
              linspace(y_range[0], y_range[1], static_cast<size_t>(mesh_density));
          auto [X, Y] = meshgrid(x, y);
          auto Z = transform(X, Y, fn);
          surface_handle l = this->mesh(X, Y, Z, {});
  
          return l;
      }
  
      /// Lambda function mesh - parametric function
      surface_handle axes_type::fmesh(fcontour_function_type funx,
                                      fcontour_function_type funy,
                                      fcontour_function_type funz,
                                      const std::array<double, 2> &u_range,
                                      const std::array<double, 2> &v_range,
                                      double mesh_density) {
          axes_silencer temp_silencer_{this};
  
          std::vector<double> u =
              linspace(u_range[0], u_range[1], static_cast<size_t>(mesh_density));
          std::vector<double> v =
              linspace(v_range[0], v_range[1], static_cast<size_t>(mesh_density));
          auto [U, V] = meshgrid(u, v);
          auto X = transform(U, V, funx);
          auto Y = transform(U, V, funy);
          auto Z = transform(U, V, funz);
          surface_handle l = this->mesh(X, Y, Z, {});
  
          return l;
      }
  
      /// Lambda function mesh
      /// Grid / Both ranges in the same array size 4
      surface_handle axes_type::fmesh(fcontour_function_type fn,
                                      const std::array<double, 4> &xy_range,
                                      double mesh_density) {
          return this->fmesh(fn, {xy_range[0], xy_range[1]},
                             {xy_range[2], xy_range[3]},
                             static_cast<size_t>(mesh_density));
      }
  
      /// Lambda function mesh - Parametric
      /// Parametric / Both ranges in the same array size 4
      surface_handle axes_type::fmesh(fcontour_function_type funx,
                                      fcontour_function_type funy,
                                      fcontour_function_type funz,
                                      const std::array<double, 4> &uv_range,
                                      double mesh_density) {
          return this->fmesh(funx, funy, funz, {uv_range[0], uv_range[1]},
                             {uv_range[2], uv_range[3]},
                             static_cast<size_t>(mesh_density));
      }
  
      /// Function mesh
      /// Grid / Both ranges in the same array size 2
      surface_handle axes_type::fmesh(fcontour_function_type fn,
                                      const std::array<double, 2> &xy_range,
                                      double mesh_density) {
          return this->fmesh(fn, xy_range, xy_range,
                             static_cast<size_t>(mesh_density));
      }
  
      /// Function mesh
      /// Parametric / Both ranges in the same array size 2
      surface_handle axes_type::fmesh(fcontour_function_type funx,
                                      fcontour_function_type funy,
                                      fcontour_function_type funz,
                                      const std::array<double, 2> &uv_range,
                                      double mesh_density) {
          return this->fmesh(funx, funy, funz, uv_range, uv_range,
                             static_cast<size_t>(mesh_density));
      }
  
      // z = f(x,y)
      using fsurf_function_type = std::function<double(double, double)>;
  
      /// Function surf - Core function
      surface_handle axes_type::fsurf(fcontour_function_type fn,
                                      const std::array<double, 2> &x_range,
                                      const std::array<double, 2> &y_range,
                                      std::string_view line_spec,
                                      double mesh_density) {
          axes_silencer temp_silencer_{this};
  
          std::vector<double> x =
              linspace(x_range[0], x_range[1], static_cast<size_t>(mesh_density));
          std::vector<double> y =
              linspace(y_range[0], y_range[1], static_cast<size_t>(mesh_density));
          auto [X, Y] = meshgrid(x, y);
          auto Z = transform(X, Y, fn);
          surface_handle l = this->surf(X, Y, Z, {}, line_spec);
  
          return l;
      }
  
      /// Function surf - Parametric
      /// Core parametric function
      surface_handle axes_type::fsurf(fcontour_function_type funx,
                                      fcontour_function_type funy,
                                      fcontour_function_type funz,
                                      const std::array<double, 2> &u_range,
                                      const std::array<double, 2> &v_range,
                                      std::string_view line_spec,
                                      double mesh_density) {
          axes_silencer temp_silencer_{this};
  
          std::vector<double> u =
              linspace(u_range[0], u_range[1], static_cast<size_t>(mesh_density));
          std::vector<double> v =
              linspace(v_range[0], v_range[1], static_cast<size_t>(mesh_density));
          auto [U, V] = meshgrid(u, v);
          auto X = transform(U, V, funx);
          auto Y = transform(U, V, funy);
          auto Z = transform(U, V, funz);
          surface_handle l = this->surf(X, Y, Z, {}, line_spec);
  
          return l;
      }
  
      /// Function surf
      /// Grid / Both ranges in the same array size 4
      surface_handle axes_type::fsurf(fcontour_function_type fn,
                                      const std::array<double, 4> &xy_range,
                                      std::string_view line_spec,
                                      double mesh_density) {
          return this->fsurf(fn, {xy_range[0], xy_range[1]},
                             {xy_range[2], xy_range[3]}, line_spec,
                             static_cast<size_t>(mesh_density));
      }
  
      /// Parametric / Both ranges in the same array size 4
      surface_handle axes_type::fsurf(fcontour_function_type funx,
                                      fcontour_function_type funy,
                                      fcontour_function_type funz,
                                      const std::array<double, 4> &uv_range,
                                      std::string_view line_spec,
                                      double mesh_density) {
          return this->fsurf(funx, funy, funz, {uv_range[0], uv_range[1]},
                             {uv_range[2], uv_range[3]}, line_spec,
                             static_cast<size_t>(mesh_density));
      }
  
      /// Function surf
      /// Grid / Both ranges in the same array size 2
      surface_handle axes_type::fsurf(fcontour_function_type fn,
                                      const std::array<double, 2> &xy_range,
                                      std::string_view line_spec,
                                      double mesh_density) {
          return this->fsurf(fn, xy_range, xy_range, line_spec,
                             static_cast<size_t>(mesh_density));
      }
  
      /// Function surf
      /// Parametric / Both ranges in the same array size 2
      surface_handle axes_type::fsurf(fcontour_function_type funx,
                                      fcontour_function_type funy,
                                      fcontour_function_type funz,
                                      const std::array<double, 2> &uv_range,
                                      std::string_view line_spec,
                                      double mesh_density) {
          return this->fsurf(funx, funy, funz, uv_range, uv_range, line_spec,
                             static_cast<size_t>(mesh_density));
      }
  
      /// Mesh - Core function
      surface_handle axes_type::mesh(const std::vector<std::vector<double>> &X,
                                     const std::vector<std::vector<double>> &Y,
                                     const std::vector<std::vector<double>> &Z,
                                     const std::vector<std::vector<double>> &C) {
          axes_silencer temp_silencer_{this};
  
          surface_handle l = std::make_shared<class surface>(this, X, Y, Z, C);
          l->palette_map_at_surface(false);
          l->hidden_3d(true);
          this->emplace_object(l);
  
          return l;
      }
  
      /// Mesh with contour
      surface_handle axes_type::meshc(const std::vector<std::vector<double>> &X,
                                      const std::vector<std::vector<double>> &Y,
                                      const std::vector<std::vector<double>> &Z,
                                      const std::vector<std::vector<double>> &C) {
          axes_silencer temp_silencer_{this};
  
          surface_handle l = std::make_shared<class surface>(this, X, Y, Z, C);
          l->palette_map_at_surface(false);
          l->contour_base(true);
          l->hidden_3d(true);
          this->emplace_object(l);
  
          return l;
      }
  
      /// Mesh with curtain
      /// Core function
      surface_handle axes_type::meshz(const std::vector<std::vector<double>> &X,
                                      const std::vector<std::vector<double>> &Y,
                                      const std::vector<std::vector<double>> &Z,
                                      const std::vector<std::vector<double>> &C) {
          axes_silencer temp_silencer_{this};
  
          surface_handle l = this->mesh(X, Y, Z, C);
          l->curtain(true);
  
          return l;
      }
  
      /// Ribbon
      surface_handle axes_type::ribbon(const std::vector<std::vector<double>> &X,
                                       const std::vector<std::vector<double>> &Y,
                                       const std::vector<std::vector<double>> &Z,
                                       const std::vector<std::vector<double>> &C,
                                       double width) {
          axes_silencer temp_silencer_{this};
  
          surface_handle l = this->surf(X, Y, Z, C);
          l->ribbons(true);
          l->ribbon_width(width);
          l->depthorder(true);
  
          return l;
      }
  
      /// Surf - Core function
      surface_handle axes_type::surf(const std::vector<std::vector<double>> &X,
                                     const std::vector<std::vector<double>> &Y,
                                     const std::vector<std::vector<double>> &Z,
                                     const std::vector<std::vector<double>> &C,
                                     std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          surface_handle l =
              std::make_shared<class surface>(this, X, Y, Z, C, line_spec);
          this->emplace_object(l);
  
          return l;
      }
  
      /// Surf with contour - Core function
      surface_handle axes_type::surfc(const std::vector<std::vector<double>> &X,
                                      const std::vector<std::vector<double>> &Y,
                                      const std::vector<std::vector<double>> &Z,
                                      const std::vector<std::vector<double>> &C) {
          axes_silencer temp_silencer_{this};
  
          surface_handle l = std::make_shared<class surface>(this, X, Y, Z, C);
          l->contour_base(true);
          l->hidden_3d(C.empty());
          this->emplace_object(l);
  
          return l;
      }
  
      /// Waterfall - Core function
      surface_handle
      axes_type::waterfall(const std::vector<std::vector<double>> &X,
                           const std::vector<std::vector<double>> &Y,
                           const std::vector<std::vector<double>> &Z,
                           const std::vector<std::vector<double>> &C) {
          axes_silencer temp_silencer_{this};
  
          surface_handle l = this->mesh(X, Y, Z, C);
          l->waterfall(true);
  
          return l;
      }
  
      /// Graph - Undirected - Core function
      network_handle
      axes_type::graph(const std::vector<std::pair<size_t, size_t>> &edges,
                       const std::vector<double> &weights, size_t n_vertices,
                       std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
  
          network_handle l = std::make_shared<class network>(
              this, edges, weights, n_vertices, line_spec);
          this->emplace_object(l);
          this->x_axis().visible(false);
          this->y_axis().visible(false);
          if (this->x_axis().limits_mode_auto()) {
              this->x_axis().limits({l->xmin(), l->xmax()});
          }
          if (this->y_axis().limits_mode_auto()) {
              this->y_axis().limits({l->ymin(), l->ymax()});
          }
  
          return l;
      }
  
      network_handle
      axes_type::graph(const std::vector<std::pair<size_t, size_t>> &edges,
                       std::string_view line_spec) {
          return this->graph(edges, vector_1d{}, 0, line_spec);
      }
  
      /// B&W image show - Core function
      matrix_handle axes_type::imshow(
          const std::vector<std::vector<unsigned char>> &gray_scale_img) {
          axes_silencer temp_silencer_{this};
  
          // create matrix in the xlim
          matrix_handle img =
              std::make_shared<class matrix>(this, gray_scale_img);
          img->always_hide_labels(true);
          this->emplace_object(img);
          this->axis(equal);
          this->color(this->parent()->color());
          this->colormap(palette::gray());
          this->box(false);
          this->grid(false);
          this->minor_grid(false);
          this->grid_front(false);
          this->color_box(false);
          this->x_axis().visible(false);
          this->y_axis().visible(false);
          this->y_axis().reverse(true);
          this->color_box_range(0, 255);
  
          // plot if not quiet
          return img;
      }
  
      /// Core RGB / RGBA imshow function
      matrix_handle axes_type::imshow(
          const std::vector<std::vector<unsigned char>> &r_channel,
          const std::vector<std::vector<unsigned char>> &g_channel,
          const std::vector<std::vector<unsigned char>> &b_channel,
          const std::vector<std::vector<unsigned char>> &a_channel) {
          axes_silencer temp_silencer_{this};
          // Create matrix in the xlim
          matrix_handle img = std::make_shared<class matrix>(
              this, r_channel, g_channel, b_channel, a_channel);
          this->emplace_object(img);
          this->axis(equal);
          this->color(this->parent()->color());
          this->colormap(palette::gray());
          this->box(false);
          this->grid(false);
          this->minor_grid(false);
          this->grid_front(false);
          this->color_box(false);
          this->x_axis().visible(false);
          this->y_axis().visible(false);
          this->y_axis().reverse(true);
          this->color_box_range(0, 255);
          return img;
      }
  
      /// Image show from 2d vectors with image channels
      matrix_handle axes_type::imshow(
          const std::vector<std::vector<std::vector<unsigned char>>> &img) {
          axes_silencer temp_silencer_{this};
          if (img.size() == 1) {
              return this->imshow(img[0]);
          } else if (img.size() == 3) {
              return this->imshow(img[0], img[1], img[2]);
          } else if (img.size() == 4) {
              return this->imshow(img[0], img[1], img[2], img[3]);
          } else {
              return matrix_handle{nullptr};
          }
      }
  
      /// Image show from filename
      matrix_handle axes_type::imshow(const std::string &filename) {
          image_channels_t image = imread(filename);
          return this->imshow(image);
      }
  
      /// Display array as image
      matrix_handle axes_type::image(const std::vector<std::vector<double>> &C,
                                     bool scaled_colorbar) {
          axes_silencer temp_silencer_{this};
  
          // create matrix in the xlim
          matrix_handle img = std::make_shared<class matrix>(this, C);
          img->always_hide_labels(true);
          this->emplace_object(img);
          this->color(this->parent()->color());
          this->box(true);
          this->color_box(false);
          if (this->children().size() == 1) {
              this->x_axis().limits({img->xmin(), img->xmax()});
              this->y_axis().limits({img->ymin(), img->ymax()});
          }
          this->y_axis().reverse(true);
          if (!scaled_colorbar) {
              this->color_box_range(0, 255);
          }
  
          // plot if not quiet
          return img;
      }
  
      /// Display 3 arrays as image
      matrix_handle
      axes_type::image(const std::vector<std::vector<double>> &r_channel,
                       const std::vector<std::vector<double>> &g_channel,
                       const std::vector<std::vector<double>> &b_channel,
                       bool scaled_colorbar) {
          axes_silencer temp_silencer_{this};
  
          // create matrix in the xlim
          matrix_handle img = std::make_shared<class matrix>(
              this, r_channel, g_channel, b_channel);
          img->always_hide_labels(true);
          this->emplace_object(img);
          this->color(this->parent()->color());
          this->box(true);
          this->color_box(false);
          if (this->children().size() == 1) {
              this->x_axis().limits({img->xmin(), img->xmax()});
              this->y_axis().limits({img->ymin(), img->ymax()});
          }
          if (!scaled_colorbar) {
              this->color_box_range(0, 255);
          }
  
          // plot if not quiet
          return img;
      }
  
      /// Show 2d-array as image
      matrix_handle axes_type::image(double x_min, double x_max, double y_min,
                                     double y_max,
                                     const std::vector<std::vector<double>> &C,
                                     bool scaled_colorbar) {
          axes_silencer temp_silencer_{this};
  
          matrix_handle img = this->image(C, scaled_colorbar);
          img->x(x_min);
          img->w(x_max - x_min + 1);
          img->y(y_min);
          img->h(y_max - y_min + 1);
          this->x_axis()
  
              .limits({img->
  
                       xmin(),
                       img
  
                           ->xmax()
  
              });
          this->y_axis()
  
              .limits({img->
  
                       ymin(),
                       img
  
                           ->ymax()
  
              });
  
          // plot if not quiet
          return img;
      }
  
      /// Show image from tuple of rgb 2d vectors
      matrix_handle
      axes_type::image(const std::tuple<std::vector<std::vector<double>>,
                                        std::vector<std::vector<double>>,
                                        std::vector<std::vector<double>>> &img) {
          return this->image(std::get<0>(img), std::get<1>(img), std::get<2>(img),
                             false);
      }
  
      /// Annotate plot with text
      /// \param x Positions
      /// \param y Positions
      /// \param texts Texts
      /// \return
      labels_handle axes_type::text(const vector_1d &x, const vector_1d &y,
                                    const std::vector<std::string> &texts) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
  
          this->next_plot_replace(false);
          labels_handle l = std::make_shared<class labels>(this, x, y, texts);
          l->alignment(labels::alignment::left);
          this->emplace_object(l);
  
          this->next_plot_replace(p2);
          return l;
      }
  
      /// Annotate plot with single text
      labels_handle axes_type::text(double x, double y, std::string_view str) {
          return this->text(std::vector{x}, std::vector{y}, std::vector{std::string(str)});
      }
  
      /// Annotate plot with same text at many positions
      labels_handle axes_type::text(const vector_1d &x, const vector_1d &y,
                                    std::string_view str) {
          return this->text(x, y, std::vector<std::string>(x.size(), std::string(str)));
      }
  
      /// Annotate plot with arrow
      vectors_handle axes_type::arrow(double x1, double y1, double x2,
                                      double y2) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
  
          double u = x2 - x1;
          double v = y2 - y1;
          vectors_handle q =
              std::make_shared<class vectors>(this, vector_1d{x1}, vector_1d{y1},
                                              vector_1d{u}, vector_1d{v}, "-k");
          q->line_width(1.5);
          this->emplace_object(q);
  
          this->next_plot_replace(p2);
          return q;
      }
  
      /// Annotate plot with single line between two points
      line_handle axes_type::line(double x1, double y1, double x2, double y2) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
  
          auto l = this->plot(std::vector{x1, x2}, std::vector{y1, y2}, "-k");
          l->line_width(1.5);
  
          this->next_plot_replace(p2);
          return l;
      }
  
      /// Annotate plot with text and arrow
      std::pair<labels_handle, vectors_handle>
      axes_type::textarrow(double x1, double y1, double x2, double y2,
                           std::string_view str) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
  
          auto t = this->text(x1, y1, str);
          double u = x2 - x1;
          if (u > 0) {
              t->alignment(labels::alignment::right);
          }
          auto q = this->arrow(x1, y1, x2, y2);
  
          this->next_plot_replace(p2);
          return std::make_pair(t, q);
      }
  
      /// Annotate plot with rectangle
      line_handle axes_type::rectangle(double x, double y, double w, double h,
                                       double curvature) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
  
          line_handle l = nullptr;
          if (curvature == 0.) {
              std::vector<double> x_ = {x, x + w, x + w, x, x};
              std::vector<double> y_ = {y, y, y + h, y + h, y};
              l = plot(x_, y_, "-k");
          } else {
              // Radius of the round section of the rectangle
              double radius = min(w, h) * 0.5 * curvature;
  
              // Rounded corners
              auto [xne, yne] = pol2cart(linspace(0, pi / 2), radius);
              auto [xnw, ynw] = pol2cart(linspace(pi / 2, pi), radius);
              auto [xsw, ysw] = pol2cart(linspace(pi, 3 * pi / 2), radius);
              auto [xse, yse] = pol2cart(linspace(3 * pi / 2, 2 * pi), radius);
  
              // Use corners
              std::vector<double> x_;
              x_ = concat(x_, transform(xne, [&](double c) {
                              return x + w - radius + c;
                          }));
              x_ = concat(
                  x_, transform(xnw, [&](double c) { return x + radius + c; }));
              x_ = concat(
                  x_, transform(xsw, [&](double c) { return x + radius + c; }));
              x_ = concat(x_, transform(xse, [&](double c) {
                              return x + w - radius + c;
                          }));
              x_.emplace_back(x_.front());
  
              std::vector<double> y_;
              y_ = concat(y_, transform(yne, [&](double c) {
                              return y + h - radius + c;
                          }));
              y_ = concat(y_, transform(ynw, [&](double c) {
                              return y + h - radius + c;
                          }));
              y_ = concat(
                  y_, transform(ysw, [&](double c) { return y + radius + c; }));
              y_ = concat(
                  y_, transform(yse, [&](double c) { return y + radius + c; }));
              y_.emplace_back(y_.front());
  
              l = plot(x_, y_, "-k");
          }
  
          this->next_plot_replace(p2);
          return l;
      }
  
      std::pair<labels_handle, line_handle>
      axes_type::textbox(double x, double y, double w, double h,
                         std::string_view str) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
  
          auto r = rectangle(x, y, w, h);
          auto t = text(x + w / 2, y + h / 2, str);
          t->alignment(labels::alignment::center);
  
          this->next_plot_replace(p2);
          return std::make_pair(t, r);
      }
  
      line_handle axes_type::fill(const vector_1d &x, const vector_1d &y,
                                  std::string_view line_spec) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
  
          auto l = this->plot(x, y, line_spec);
          l->fill(true);
  
          this->next_plot_replace(p2);
          return l;
      }
  
      line_handle axes_type::ellipse(double x, double y, double w, double h) {
          axes_silencer temp_silencer_{this};
          bool p2 = this->next_plot_replace();
          this->next_plot_replace(false);
  
          auto t = linspace(0, 2 * pi);
          auto [x_, y_] = elliptic2cart(t, w / 2., h / 2.);
          x_ = transform(x_, [&](double x_) { return x_ + w / 2 + x; });
          y_ = transform(y_, [&](double y_) { return y_ + h / 2 + y; });
          auto l = plot(x_, y_, "-k");
  
          this->next_plot_replace(p2);
          return l;
      }
  
      /// Create boxplot from vector 1d
      box_chart_handle axes_type::boxplot(const std::vector<double> &data) {
          box_chart_handle l = std::make_shared<class box_chart>(this, data);
          this->emplace_object(l);
          return l;
      }
  
      /// Create boxplot from vector 2d
      box_chart_handle
      axes_type::boxplot(const std::vector<std::vector<double>> &data) {
          std::vector<double> y_data;
          std::vector<double> x_data;
          size_t index = 1;
          for (auto it_i = data.begin(); it_i != data.end(); ++it_i) {
              for (auto it_j = it_i->begin(); it_j != it_i->end(); ++it_j) {
                  y_data.emplace_back(*it_j);
                  x_data.emplace_back(static_cast<double>(index));
              }
              ++index;
          }
          box_chart_handle l =
              std::make_shared<class box_chart>(this, y_data, x_data);
          this->emplace_object(l);
          this->x_axis().tick_values(iota(1, 1, data.size()));
          return l;
      }
  
      /// Create boxplot with groups
      box_chart_handle axes_type::boxplot(const std::vector<double> &data,
                                          const std::vector<double> &groups) {
          box_chart_handle l =
              std::make_shared<class box_chart>(this, data, groups);
          this->emplace_object(l);
          return l;
      }
  
      /// Create boxplot with group strings
      box_chart_handle
      axes_type::boxplot(const std::vector<double> &y_data,
                         const std::vector<std::string> &groups) {
          // map<group,index>
          std::map<std::string, double> category_indexes;
          std::vector<double> group_indexes;
          for (const std::string &group : groups) {
              // look for group
              auto it = category_indexes.find(group);
              // if already there
              if (it != category_indexes.end()) {
                  // keep its numeric equivalent
                  group_indexes.emplace_back(it->second);
              } else {
                  // put in the set and get number
                  double group_num = category_indexes.size() + 1;
                  category_indexes[group] = group_num;
                  group_indexes.emplace_back(group_num);
              }
          }
          box_chart_handle h = this->boxplot(y_data, group_indexes);
          this->xticks(iota(1, category_indexes.size()));
  
          std::vector<std::string> unique_categories(category_indexes.size());
          for (const auto &[category, index] : category_indexes) {
              unique_categories[static_cast<size_t>(round(index) - 1)] = category;
          }
          this->x_axis().ticklabels(unique_categories);
          this->x_axis().limits(
              {0.5, static_cast<double>(category_indexes.size()) + 0.5});
          return h;
      }
  
      void axes_type::draw_path(const std::vector<double> &x,
                                const std::vector<double> &y,
                                const std::array<float, 4> &color) {
          // we still have to make limits calculate and return the
          // automatic limits rather than the default limits
          auto xlimits = xlim();
          auto ylimits = ylim();
          // clamp
          std::vector<double> cx = transform(
              x, [&](double x) { return std::clamp(x, xlimits[0], xlimits[1]); });
          std::vector<double> cy = transform(
              y, [&](double y) { return std::clamp(y, ylimits[0], ylimits[1]); });
          // normalize
          auto [w,h,lm,rm,bm,tm] = calculate_margins();
          double view_width = parent_->backend_->width();
          double view_xmin = lm * view_width;
          double view_xmax = rm * view_width;
          double view_height = parent_->backend_->height();
          double view_ymin = bm * view_height;
          double view_ymax = tm * view_height;
          auto xrange = xlimits[1] - xlimits[0];
          auto yrange = ylimits[1] - ylimits[0];
          for (auto &v : cx) {
              v -= xlimits[0];
              v /= xrange;
              v *= view_xmax - view_xmin;
              v += view_xmin;
          }
          for (auto &v : cy) {
              v -= ylimits[0];
              v /= yrange;
              v *= view_ymax - view_ymin;
              v += view_ymin;
          }
          // draw the normalized path to the backend
          parent_->backend_->draw_path(cx,cy,color);
      }
  
  } // namespace matplot