Guest Post On Fluent C++

On September 11th, 2020, I had a guest post published on Jonathan Bocarra’s blog, Fluent C++. This is one of the C++ blogs I read regularly so I was very excited to have something published there.

The post is titled “Replacing CRTP Static Polymorphism With Concepts” and in it I iteratively update legacy CRTP code with C++20 concepts. This is important because you may not be able to switch your whole code base to a concept based design all at once. Even so, you may find that a partial refactor will allow you to take advantage of the expressiveness and improved compile time error messages provided by concepts.

You can check out the post here.

Managing Class Options

I’ve found myself recently needing to handle passing and managing options across all of the components of our application.
I wanted all components to have their own options that could be managed externally.

Approach

I nest the options struct in the class and let it be qualified by the class name.

class MyComponent {
public:
  struct Options {
    struct {
      bool some_option{false};
      /* some other options */
    } MyComponent;
  };
  
  MyComponent() {}
  MyComponent(std::shared_ptr<Options>& options) :
      m_options{options} {}
  /* ... */
private:
  std::shared_ptr<Options> m_options{new Options};
};

We will devise a way to aggregate options and get a specific slice of options from it.

template <typename... TOptions>
struct OptionsAggregate : public TOptions... {
  using Type = OptionsAggregate<TOptions...>;
  template <typename TOptionType>
  auto Get() -> TOptionType& {
    static_assert(std::is_base_of_v<TOptionType, Type>,
      "Requested options are not a base of this aggregate");
    return static_cast<TOptionType&>(*this);
  }
};

So this variadic class template will take several options classes and inherit from them. It also has a Get method to get a slice of the aggregate.

We will have a ComponentManager that manages all of the components and their options.

template <typename TOptions>
class ComponentManager {
public:
  ComponentManager() {}
  ComponentManager(std::shared_ptr<MyComponent> component,
                   std::shared_ptr<Options> options) :
     m_options{std::move(options)} {}
  void Run();
private:
  std::shared_ptr<TOptions> m_options{new TOptions};
};

We can then use all of these things like the following:

struct GeneralOptions {
  bool debug{false};
};

using ComponentOptions = OptionsAggregate<GeneralOptions, MyComponent::Options>;
auto options = std::make_shared<ComponentOptions>();

options->Get<GeneralOptions>().debug = true;
/* set other options */

// the options that component can see are only a slice of the overall options
auto component = std::make_shared<MyComponent>(options);

ComponentManager manager{component, options};

A problem I discovered with this method appears when using the slicing syntax from ComponentManager methods which results in very unfortunate syntax. Imagine component manager had a method called Run:

template<typename TOptions>
void ComponentManager<TOptions>::Run() {
  auto& general_options = m_options->template Get<GeneralOptions>();
  /* use general_options */
}

This occurs because m_options depends on a template parameter as well as OptionsAggregate::Get.

Check Compiler Explorer for a full example.

Generic Visitor Builder

Recently, I’ve been implementing an internal scene graph data structure and I want to show my solution for building generic scene graph visitors.

Lets consider a basis type for our tree.

class base_node {
public:
  base_node() : m_name{"default"} {}
  base_node(std::string name) : m_name{std::move(name)} {}

  [[nodiscard]] auto& name() { return m_name; }
  [[nodiscard]] auto& name() const { return m_name; }
  [[nodiscard]] auto& children() { return m_children; }
  [[nodiscard]] auto& children() const { return m_children; }

private:
  std::string m_name;
  std::vector<base_node> m_children;
};

We can create a new scene graph like this:

  auto layer_1 = base_node{"parent"};
  
  auto& layer_2 = layer_1.children();
  layer_2 = {{"child_1"},{"child_2"},{"child_3"},{"child_4"},{"child_5"}};

  auto& layer_3 = layer_2[1].children();
  layer_3 = {{"child_6"},{"child_7"},{"child_8"},{"child_9"},{"child_10"}};

It has a name and a vector of base_nodes to hold its children.

So now lets say you want to apply some function to the whole scene graph; lets say you want to print the hierarchy.

We could write a visitor functor which takes a base_node and recursively calls itself with all of its children printing the name of each node. For this example, we will visit in a preorder fashion.

We will also add some state to make it look pretty.

struct print_hierarchy {
  void operator()(base_node& parent) {
    auto put_indent = [count = m_level] (const char indent_char) mutable {
      while(count--) std::cout << indent_char;
    };
    auto increase_indent = [&] () {
      m_level += 2;
    };
    auto decrease_indent = [&] () {
      m_level -= 2;
    };

    put_indent('-');

    std::cout << parent.name() << '\n';

    increase_indent();
    for(auto& node: parent.children()) {
      (*this)(node);
    }
    decrease_indent();
  }
  
private:
  size_t m_level{0};
};

If we call this on our hierarchy, we get:

parent
--child_1
--child_2
----child_6
----child_7
----child_8
----child_9
----child_10
--child_3
--child_4
--child_5

We can make this more generic by replacing the cout with a call to some callable.

We need to pass a Callable template argument to our visitor. We will also rename it to hierarchy_visitor

struct hierarchy_visitor {
  template<class TCallable>
  void operator()(base_node& parent, TCallable callable) {
    auto put_indent = [count = m_level] (const char indent_char) mutable {
      while(count--) std::cout << indent_char;
    };
    auto increase_indent = [&] () {
      m_level += 2;
    };
    auto decrease_indent = [&] () {
      m_level -= 2;
    };

    put_indent('-');

    callable(parent);

    increase_indent();
    for(auto& node: parent.children()) {
      (*this)(node, callable);
    }
    decrease_indent();
  }
private:
  size_t m_level{0};
};

We can now do this:

auto print_node = [](base_node& parent) {
  std::cout << parent.name() << '\n';
};

auto visitor = hierarchy_visitor();
visitor(layer_1, print_node);

But now we have a new problem. We have formatting hard coded into our visitor. It isn’t as generic as it could be.

We want to refactor out the formatting code. We can move put_indent to the callable and make the callable accept a layer parameter that tells the current level in the hierarchy visitor. It makes since for a hierarchy visitor to know its current level so level will remain inside the visitor. We will also rename the lambdas to something that sounds more generic.

struct hierarchy_visitor {
  template <class TCallable>
  void operator()(base_node& parent, TCallable callable) {
    auto increase_level = [&]() { m_level += 1; };
    auto decrease_level = [&]() { m_level -= 1; };

    callable(parent, m_level);

    increase_level();
    for (auto& node : parent.children()) {
      (*this)(node, callable);
    }
    decrease_level();
  }

 private:
  size_t m_level{0};
};

auto print_node = [indent_char = '-', indent_width = 2](base_node& parent, auto level) {
  auto put_indent = [&indent_width, level](auto indent_char) mutable {
    level *= indent_width;
    while (level--) std::cout << indent_char;
  };

  put_indent(indent_char);
  std::cout << parent.name() << '\n';
};

What if we don’t want to pass our callable in every recursive call? We can allow our hierarchy visitor to store it as a private member.

template<class TCallable>
struct hierarchy_visitor {
  void operator()(base_node& parent) {
    auto increase_level = [&]() { m_level += 1; };
    auto decrease_level = [&]() { m_level -= 1; };

    m_callable(parent, m_level);

    increase_level();
    for (auto& node : parent.children()) {
      (*this)(node, callable);
    }
    decrease_level();
  }

 private:
  TCallable m_callable;
  size_t m_level{0};
};

But now we have an issue where we cant determine the type of a lambda to instantiate an instance of hierarchy_visitor. We will have to give hierarchy_visitor a constructor which takes a TCallable.

But even this is not enough. Since
We will have to create a builder.

template <typename Callable>
hierarchy_visitor<Callable> build_hierarchy_visitor(Callable callable) {
  return {callable};
}

Now we can build a hierarchy_visitor but what about other types of visitors?

We can use some lambda magic with variable templates to accept a generic TCallable and TVisitor.

template <template <class> class TVisitor>
auto build_visitor =
    [](auto callable) -> TVisitor<std::decay_t<decltype(callable)>> {
  return {callable};
};

Unfortunately this doesn’t work in MSVC without -std:c++latest so we’ll have to try something else for the most portable code.

template <template <class> class TVisitor, class TCallable>
TVisitor<TCallable> build_visitor(TCallable callable) {
  return {callable};
}

We can then create new visitors and call them like this:

  auto visitor = build_visitor<hierarchy_visitor>(print_node);
  visitor(layer_1);

I tested this one on Clang 9, GCC 9.2, and MSVC 19.24 and it works fine.

Here is an alternative which may be more flexible and also works on all platforms:

template <template <class> class TVisitor>
struct build_visitor_impl {
  template <typename TCallable>
  constexpr auto operator()(TCallable callable) const
      -> TVisitor<decltype(callable)> {
    return {callable};
  }
};

template <template <class> class TVisitor>
inline constexpr auto build_visitor = build_visitor_impl<TVisitor>{};

View the full code on Compiler Explorer.

Thanks for reading!

Solving Coding Challenges - UVA 823 - Number 1 Runtime

This programming challenge is called Stopper Stumper. This, as I see it, is a computational geometry problem. It is a variation of the packing problem and the simplified problem statement can be stated as such:

Determine whether three circles of random size can fit within a random triangle.

Of course, it cannot be quite so simple as this. “Stoppers” have two sides with different diameters. The problem becomes:

Determine if, in any orientation, three stoppers with random top and bottom diameters can 
fit within a random triangle.

This addition, however, is a red herring. This doesn’t really make the problem much harder. This only means that you must check both “sides” to determine the fit.

Input

From the problem description:

"The input consists of a sequence of triangle specifications and descriptions of three   
stoppers for each triangle. Each triangular space is specified by three positive integers
representing the lengths of the three sides; only valid triangles will appear in the input. A
pair of positive real numbers represents each stopper. The first number in the pair represents 
the diameter of the smaller cylinder, and the second represents the diameter of the larger 
cylinder. The final line of the input file contains zeros for all the data values."

An example:

10 10 10 2.0 3.0 1.0 2.0 1.5 3.5

Output

Triangle number 2:
All three stoppers will fit in the triangular space

How I solved it

Some observations about the problem in general:

1. Triangles are given in side lengths
2. The triangle lengths are given integers
3. There are always three stoppers
4. Stoppers have two sides
5. The stopper diameters are given as floats

So I need a few things in order to store my input data:

1. Something to hold a stopper
2. Something to hold a triangle

I started by creating a type to hold a stopper and a fixed size std::array to hold 3 stoppers. I also defined some constants that I will need later:

const float PI = 3.14159265358979;
const float EPSILON = 1E-9;

typedef std::array<float, 2> stopper; // stoppers have two sides
static std::array<stopper, 3> stopper_list;  // there are always three stoppers

I decided to use scanf instead of cin to read the input because it is reasonably fast, but I could also just turn off sync_with_stdio for faster cin. I put this in a while loop and we want to read until we encounter a triangle with 0 length sides. Even though the triangle sides are integers, I read them as floats so I don’t have to do a cast later. This indicates our end condition. I also define a count to use when printing the output:

int count = 1;
while (scanf(" %f %f %f %f %f %f %f %f %f ",
             &triangle[0], &triangle[1], &triangle[2],
             &stopper_list[0][0], &stopper_list[0][1],
             &stopper_list[1][0], &stopper_list[1][1],
             &stopper_list[2][0], &stopper_list[2][1])
         && triangle[0] != 0) { }

I sort the side lengths because it is important that the orientation of the triangle is consistent for every case:

std::sort(triangle.begin(), triangle.end());

I determined that it wasn’t very useful to have the triangles represented by side lengths so I decided to determine the angles of the triangles. I created another fixed size array to hold the angles in degrees as floats:

static std::array<float, 3> angles;

There are two common cases here for determining the angles: the normal case that involves lots of trig, and the “easy” case where the triangle is equilateral. Equilateral triangles’ angles are all 60 degrees. This is important because floating point trig operations are expensive:

if (triangle[0] == triangle[1] && triangle[1] == triangle[2]) {
  angles[0] = angles[1] = angles[2] = 60.0;  // equilateral triangle, angles always 60
}

It had been a while since I had done trig so I went to google to determine how to get the angles of a triangle from the side lengths alone.
The Law of Cosines allows us to find angles in a triangle only using side lengths: Law of Cosines

c2 = a2 + b2 − 2ab cos(C)

That’s great! So I put that in code and now I have one angle from the triangle in radians. So now we have:

if (triangle[0] == triangle[1] && triangle[1] == triangle[2]) {
  ...
} else {
  angles[0] = acos(
      (triangle[0] * triangle[0] + 
       triangle[1] * triangle[1] - 
       triangle[2] * triangle[2] / 
       (2.0 * triangle[0] * triangle[1]));  // law of cosines for first angle
}

Instead of using the Law of Cosines for every angle, I switched to the Law of Sines. We can determine an angle with two side lengths and a single angle. It also requires less operations to determine the angle and we can use it now that we have one of the angles. Law of Sines

a / sin A = b / sin B = c / sin C

Now we have a second angle. We actually want degrees so I convert the angles to degrees:

...
} else {
  angles[0] = acos(
      (triangle[0] * triangle[0] + 
       triangle[1] * triangle[1] - 
       triangle[2] * triangle[2] / 
       (2.0 * triangle[0] * triangle[1]));  // law of cosines for first angle
       
  angles[1] = asin(triangle[0] / triangle[2] * sin(angles[0]))
                  * (180.0 / PI);           // law of sines for second angle

  angles[0] *= (180.0 / PI);                // convert angle 0 to degrees
}

The last angle is the easiest since the degrees of a triangle always adds up to 180 degrees. The angles are sorted so they correspond to the side length opposite of them.:

...
   angles[2] = 180.0 - (angles[0] + angles[1]);  // subtraction for last angle
   std::sort(angles.begin(), angles.end());
...

I now wanted to map the triangle to the Cartesian plane. I created a type to hold grid points and a fixed size array to hold 3 grid_points:

struct grid_point {
  grid_point(float x = 0.0, float y = 0.0) : x(x), y(y) {}
  float x, y;
};

static std::array<grid_point, 3> vertices;

Now that we have the side lengths and the angles for this triangle, we can map it to the Cartesian plane.
I start by always mapping the vertex with the smallest angle to (0,0):

vertices[0].x = 0.0;            // vertices for 0 points always 0
vertices[0].y = 0.0;            // vertices for 0 points always 0

And the vertex with the largest angle onto the x-axis:

vertices[1].x = triangle[2];  // vertices for second point made up of (triangle[2], 0.0)
vertices[1].y = 0.0;

The third point is a little trickier than the others since they are on the axis. In order to determine the coordinates for the last point I used the following formulas represented in code:

vertices[2].x = cos(angles[0] * PI / 180.0) * triangle[1];
vertices[2].y = sin(angles[0] * PI / 180.0) * triangle[1];

We now have all of the information we need for the triangle. It is time to move on to the stoppers. We have to determine the optimal stopper to “flip” before we can check the fit. The algorithm I created to determine the best stopper to flip is as follows:

If the large side of all three stoppers is the same size
  Flip the stopper with the minimal small side
Else
  Flip the stopper with the maximal large side

In code this looks like the following:

size_t index;
if (stopper_list[0][1] == stopper_list[1][1] && 
    stopper_list[1][1] == stopper_list[2][1]) {
  index = std::distance(stopper_list.begin(), std::min_element(
            stopper_list.begin(), stopper_list.end(), 
            [] (const stopper& a, const stopper& b) { 
              return a[0] < b[0];
            }));
} else {
  index = std::distance(stopper_list.begin(), std::max_element(
            stopper_list.begin(), stopper_list.end(),
            [](const stopper& a, const stopper& b) { 
              return a[1] < b[1];
            }));
}
std::swap(stopper_list[index], stopper_list[0]); // the stopper in position 0 is 'flipped'

Now we want to go through the 6 permutations of stopper positions. I hard code a const array of stopper positions to loop through.

static const std::array<std::array<int, 3>, 6> possibilities{
  0, 1, 2,
  0, 2, 1,
  1, 0, 2,
  1, 2, 0,
  2, 0, 1,
  2, 1, 0 };

An easy short circuit case to check is to see if all of the circles can fit inside the triangle in any orientation by themselves. The maximum radius circle that could possibly fit inside the triangle is given by the formula:

max_radius = (2 * triangle_area) /  triangle_perimeter

I calculate the triangle area with the Heron’s formula Heron’s Formula:

A = sqrt( s ( s-a ) ( s-b ) ( s-c ) )

Where s is the semiperimeter and a, b, and c are the side lengths of the triangle. The semiperimeter is the perimeter of the triangle divided by 2. In code this looks like this:

double tri_perim = (triangle[0] + triangle[1] + triangle[2]);
double tri_semiperim = tri_perim / 2.0;
double tri_area = sqrt(tri_semiperim * 
                    (tri_semiperim - triangle[0]) * 
                    (tri_semiperim - triangle[1]) * 
                    (tri_semiperim - triangle[2]));

The maximum radius that can fit in the triangle is given by:

double max_rad = (2.0 * tri_area) / tri_perim;

I then define a bool called possible with a default value of true and loop over the stoppers large side to check for ones that can’t possibly fit:

bool possible = true;
for (int j = 0; j < 3; j++) {
  if ((stopper_list[j][1] / 2.0) > max_rad) {
    possible = false;
    break;
  }
}

I then check if possible is true and if it isn’t, I print out the negative answer. I define a bool named good with a default value of false to use in our true case:

bool good = false;
if (possible) {
  ...
}
printf("Triangle number %d:\n", count);
count++;
if (good) {
  puts("All three stoppers will fit in the triangular space");
} else {
  puts("Stoppers will not fit in the triangular space");
}

Now we have to implement the fit for the stoppers in the true case. We want to loop through the 6 possible permutations and check the placement of each circle. I define three grid_point object to use as center points for our stoppers.

for (int i = 0; i < 6; i++) {
  grid_point a, b, c;
}

I then put the circles with respect to the current permutation in the corners in order. For the first grid point a, we put the circle in the (0,0) corner of the triangle. In order to get the center for this circle I use this formula:

a_x = radius_first_circle / tan (THETA_1)
a_y = radius_first_circle

```c++
a.x = (stopper_list[possibilities[i][0]][1] / 2.0) /
      tan((angles[0] / 2.0) * PI / 180.0);          // get x of center point of first circle
a.y = (stopper_list[possibilities[i][0]][1] / 2.0); // get y of center point of first circle

The subscript operators get pretty hard to read here. This says get the diameter of the large side of the circle in position 0 in the current permutation and divide it by 2.0 to get the radius, then divide it by the tangent of the smallest angle at position 0 in our angles array.
For the second circle at corner ( triangle[2] (the longest side) , 0.0 ), we use a similar formula to the point a:

b_x = longest_side - (radius_second_circle / tan (THETA_2)
b_y = radius_second_circle

Here is the same thing in code:

b.x = triangle[2] - ((stopper_list[possibilities[i][1]][1] / 2.0) /
                    tan((angles[1] / 2.0) * PI / 180.0));   // get x of center point of second circle
b.y = (stopper_list[possibilities[i][1]][1] / 2.0);         // get y of center point of second circle

Calculating the last circle center is the worst of all. The approach I took to calculate it was to calculate x and y as if they were at (0,0) and then rotate it to the correct angle and add the x, y offset of the corner point. The initial cx and cy calculation are the same as calculating a:

tmp_cx = radius_third_circle / tan(THETA_3)
tmp_cy = radius_third_circle

And in code:

float tmp_cx = (stopper_list[possibilities[i][2]][1] / 2.0) /
               tan((angles[2] / 2.0) * PI / 180.0);           // get x of center point of first circle
float tmp_cy = (stopper_list[possibilities[i][2]][1] / 2.0);  // get y of center point of first circle

At this point we need to calculate a rotation angle. The rotation angle is 180 - ( 90 - the angle of the left side of the triangle as if it were a right triangle). This is the formula:

rot = 180.0 - (90.0 - (acos(y_component_last_point / median_side) * 180.0 / PI))

And the code:

float rotation_angle = 180.0 + (90.0 - (acos(vertices[2].y / triangle[1]) * 180.0 / PI));

Then we perform the rotation and add the offset of the corner point to get the final coordinates for c:

c.x = ((tmp_cx)*cos((rotation_angle)* PI / 180.0))
      - ((tmp_cy)*sin((rotation_angle)* PI / 180.0)) + vertices[2].x;
c.y = ((tmp_cy)*cos((rotation_angle)* PI / 180.0))
      + ((tmp_cx)*sin((rotation_angle)* PI / 180.0)) + vertices[2].y;

I then use an algorithm do do a few things:

if a is the 'flipped' stopper
  tmp_index_1 = 1
if b is the 'flipped' stopper
  tmp_index_2 = 1

tmp_lhs = distance between a and b
tmp_rhs = radius(a) + radius(b)

if (tmp_lhs <= (tmp_rhs + EPSILON))
  a and b collide

For all pairs I check if one of them is the flipped stopper and change indexes accordingly. I then get the distance between the stopper centers using the distance formula. The distance formula is Distance Formula:

c = sqrt( ( x_a - x_b ) ^2 + ( y_a - y_b ) ^2 ))

These blocks are as follows:

tmp_index_1 = tmp_index_2 = 1;
if (possibilities[i][0] == 0) {
  tmp_index_1 = 0;
} else if (possibilities[i][1] == 0) {
  tmp_index_2 = 0;
}
tmp_lhs = sqrt(pow(b.x - a.x, 2.0) + pow(b.y - a.y, 2.0));
tmp_rhs = stopper_list[possibilities[i][0]][tmp_index_1] / 2.0
          + stopper_list[possibilities[i][1]][tmp_index_2] / 2.0;
if (tmp_lhs <= (tmp_rhs + EPSILON)) {
  continue;
}

I check stoppers 0 vs 1, 1 vs 2, and 2 vs 0 for the front case. I then do the same for the back case. If all 6 ases are passed, good is set to true and the loop short circuits.

If after the loop, good is true, the positive case is printed. This code was able to pass all test cases and get the first 0.000s runtime on UVA on 2019-02-01 Statistics.

I also implemented the same code in Python and used PIL to render the cases. gist

I really enjoy these computational geometry problems. You get a good opportunity to use many mathematical techniques and I hope to do more like this in the future.

View the full source code on github