Thursday, September 12, 2019

Solving Coding Challenges - UVA 823 - Number 1 Runtime

UVA 823 Blog
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

No comments:

Post a Comment

Note: Only a member of this blog may post a comment.