Showing posts with label distance. Show all posts
Showing posts with label distance. Show all posts

Wednesday, November 10, 2021

C++20 concepts and traits

Intro

I like writing libraries. In 2007 I refactored the heart of our library using types into a library using templates. This looked a bit like:

struct point { double x, y; };

and then some functions using these types (there was also a polygon, etc).

When I was satisfied with the results I submitted this as a preview to Boost. The general reaction was: interesting, but it is not concept based, and therefore not generic enough.

So I had to make the types concept based, but I didn't know what it was. And I found out that concepts were not yet supported by the language - so how could I use them... That became all clear, with help of reviewers and joiners. The library was refactored again into a concept based library, it was accepted and it is now known as Boost.Geometry.

Now (14 years later) the C++ language has support for concepts.

This blog is targeted to the approach Boost.Geometry uses. It uses traits to adapt normal structs or classes to concepts. Functions in the library use those concepts. 

The example in  this blog

We start again with a point, and a function calculating the distance:

struct point { double x, y; };
template <typename Point>
double distance_using_point(const Point& p1, const Point& p2)
{
  auto sqr = [](const double v) -> double { return v * v; };
  return std::sqrt(sqr(p1.x - p2.x) + sqr(p1.y - p2.y));
}


This function is template based but not generic. It only works for any point type having .x and .y member variables. But suppose it is not a member variable but a function .x(). Or you are using a std::array with two elements [0] and [1] (denoting x, y).

Traits as free functions

We will use traits for that. There are more ways to define and use traits. A first, concise way (not recommended! see below for a better version) is to just define free templated functions:
 
namespace point_traits
{
  template <typename P> double get_x(const P& p);
  template <typename P> double get_y(const P& p);
}

and these functions can then be specialized for the type:
 
namespace point_traits
{
  template<> inline double get_x<point>(const point& p) { return p.x; }
  template<> inline double get_y<point>(const point& p) { return p.y; }
}

And this works! Suppose you are to write a distance function, using these traits, it could look like this:

template <typename P1, typename P2>
double
calculate_distance(const P1& p1, const P2& p2)
{
  auto sqr = [](const double x) -> double { return x * x; };
  return std::sqrt(sqr(point_traits::get_x<P1>(p1) - point_traits::get_x<P2>(p2))
                 + sqr(point_traits::get_y<P1>(p1) - point_traits::get_y<P2>(p2)));
}

If you use a std::array, then you can specialize it too:

namespace point_traits
{
  template<> inline double get_x<std::array<double, 2>>(const std::array<double, 2>& p) { return p[0]; }
  template<> inline double get_y<std::array<double, 2>>(const std::array<double, 2>& p) { return p[1]; }
}
 

And you can calculate the distance of two different types, for example:

  point mp{1, 1};
  std::array<double, 2> sa{2, 2};
  std::cout << calculate_distance(mp, sa) << std::endl;


It writes: 1.41421 as expected.

 

 

Evaluation of traits-as-free-functions 


These traits give the example already a kind of concept based look. And it works, but there are some problems:

  • functions cannot be specialized partially - this can be inconvenient
  • suppose we want floats too, and long doubles, and we want to have that coordinate type defined in a meta-function. Where to define it?
  • suppose we forget a traits adaptation: there is not appearing any compiler error message! Because it knows the generic templated function, and will not complain. It gives a linker error.
  • it is (therefore) not suitable for C++20 concepts 

The linker error, in case the specialization for point is forgotten or not right, will look like:

undefined reference to `double point_traits::get_x<point>(point const&) 
undefined reference to `double point_traits::get_y<point>(point const&) 
 
The linker just can't find these definitions. Compilation was fine. You can live with this, but there is a better alternative.

Concepts 

Now we will try to make a C++ 20 concept using these traits:

template <typename P> 
concept IsPoint = requires (P p) 
{ 
  point_traits::get_x<P>(p);
  point_traits::get_y<P>(p);
};

template <typename P1, typename P2>
requires IsPoint<P1> && IsPoint<P2>
double calculate_distance(const P1& p1, const P2& p2)
{
  // the same implementation
}


So we added, in bright magenta:
  1. a new concept IsPoint. That concept tries to use the functions in namespace point_traits. If they are not there, the type P does not fulfill the concept.
  2. a require clause decorating the distance function, stating that both point types should fulfill the IsPoint concept.

How cool is that!

However, even using this concept, there is the same linker error message. There is no neat compiler message about concepts, which we wished and maybe expected...

 

Traits as a struct 

To fix this, the traits can be modeled as structures. So we adapt the code, we introduce a structure (here replacing the namespace - but that is only for the sake of the sample - these struct traits are often placed in an own namespace). And the structure is short!

template <typename P> struct point_traits {};

So it is an empty structure. The magic happens in the specializations:

template<> struct point_traits<point> 
{
  static double get_x(const point& p) { return p.x; }
  static double get_y(const point& p) { return p.y; }
};

template<> struct point_traits<std::array<double, 2>> 
{
  static double get_x(const std::array<double, 2>& p) { return p[0]; }
  static double get_y(const std::array<double, 2>& p) { return p[1]; }
};

It looks nearly the same, just that now the struct is templatized with point or the std::array, instead of the free functions. And the free functions are now static member functions.

The concept is also similar, it's just the <P> part moving a bit left to the structure:

template <typename P> 
concept IsPoint = requires (P p) 
{ 
  point_traits<P>::get_x(p); 
  point_traits<P>::get_y(p); 
};

And so does the final generic concept-based distance function:

template <typename P1, typename P2>
requires IsPoint<P1> && IsPoint<P2>
double calculate_distance(const P1& p1, const P2& p2)
{
  auto sqr = [](const double x) -> double { return x * x; };
  return std::sqrt(sqr(point_traits<P1>::get_x(p1) - point_traits<P2>::get_x(p2)) 
                 + sqr(point_traits<P1>::get_y(p1) - point_traits<P2>::get_y(p2)));
}

And if we now forget a specialization? We get a neat error message. In clang (13) it reads:

prog.cc:62:16: error: no matching function for call to 'calculate_distance'
  std::cout << calculate_distance(mp, sa) << std::endl;
               ^~~~~~~~~~~~~~~~~~
prog.cc:45:8: note: candidate template ignored: constraints not satisfied [with P1 = point, P2 = std::array<double, 2>]
double calculate_distance(const P1& p1, const P2& p2)
       ^
prog.cc:43:10: note: because 'point' does not satisfy 'IsPoint'
requires IsPoint<P1> && IsPoint<P2>
         ^
prog.cc:35:20: note: because 'point_traits<P>::get_x(p)' would be invalid: no member named 'get_x' in 'point_traits<point>'
  point_traits<P>::get_x(p); 
                   ^
1 error generated. 
 
You can use C++20 and concepts, for example, in Wandbox.

Samples from libraries

Boost.Geometry uses this approach, also generalized for coordinate type and dimension. Here is the adaptation of the get function for std::array (in namespace boost::geometry::traits). It's partially specialized.
 
template <typename CoordinateType, std::size_t DimensionCount, std::size_t Dimension>
struct access<std::array<CoordinateType, DimensionCount>, Dimension>
{
    static inline CoordinateType get(std::array<CoordinateType, DimensionCount> const& a)
    {
        return a[Dimension];
    }
};

MapBox uses a similar approach, here is the adaptation for our point (in namespace mapbox::util):

template <> struct nth<0, point> 
{
  inline static auto get(const point &t) { return t.x; }
};
template <> struct nth<1, point> 
{
  inline static auto get(const point &t) { return t.y; }
};

Summary

The blog above shows, shortly, how to define traits, preferably as structs, and how to define a concept and use this concept in a free function distance.

Monday, July 26, 2010

C# Specializations by Traits IV


C# Specializations by Traits IV


Geometry

In last three blogs (I'm productive these days in blogs, instead of klocs) I explained how to to specializations in C#. This blog only explains how it begon, give different possibilities to design things, and wraps it up.


Towards specializations

Developing Boost.Geometry, in C++, and doing some more C# this year as well, I wondered how a generic geometry library in C# would look like. Ideally it would be as Boost.Geometry, of course... But Boost.Geometry is based on traits and C++ specializations, partial specializations, tag dispatching. All modern C++ generic programming stuff. Not easily doable in C#.

Calculating the distance between two arbitrary point types would already be something. I tried many methods, and discussed them with my collegue Paul den Dulk, and he helped me, suggesting more methods and some essential pieces. Putting it together, adding the dictionary as the runtime register instead of the compiletime C++ specializations, the point distance function is looking exactly the same as Boost.Geometry now!

Some notes

There are still many differences from this system with C++ specializations, of course.
  • in C++ it is compile-time and checked at compile-time
  • therefore in C++ it is very fast
  • in C# the dictionary is runtime and checked at runtime
  • because the dictionary lookup and necesary cast it is slower in C#
  • C++ has partial specializations, don't know (yet) how to do that in C#
  • C++ has type calculations, the whole template metaprogramming is based on that. 
  • in C# it is only the methods
  • because in C# it is runtime, it gives some advantages as well: you can register different traits classes when needed
  • and so, e.g. the ORM can handle the same model in different ways

Timings

As said the dictionary consultation, together with the cast, take some time. Using the geometry case distance, I timed it and it is certainly measurable.

Ten million distance calculations on my computer:
  • Using interfaces: 0.10 seconds
  • Using delegates: 0.15 seconds
  • Using traits: 0.13 seconds
  • Using a register: 0.42 seconds
So, in fact 4 times slower. It is much in this case, but this is compared to a small calculation (Pythagoras). If applied using GUI, database, disk io, or occasionally (and not 10 million times), the lookup will probably fall away.

Complete program

The program to research, with 10 possible implementations, with interfaces, wrappers, delegates, traits, registers and lambdas, is shown below.

using System;
using System.Diagnostics;

using System.Collections.Generic;


namespace GenericGeometry
{
    interface IPoint
    {
        double x();
        double y();
        void increase();
    }

    delegate double get_x<Point>(Point p);
    delegate double get_y<Point>(Point p);

    public interface ITraitsBase {}

    public interface ITraits<Point> : ITraitsBase
    {
        double get_x(Point p);
        double get_y(Point p);
    }

    class TraitsRegister
    {
        public static IDictionary<Type, ITraitsBase> register;

        public static void Register<PointType>(ITraits<PointType> t)
        {
            if (register == null)
            {
                register = new Dictionary<Type, ITraitsBase>();
            }
            register[typeof(PointType)] = t;
        }
    }

    class LambdaRegister
    {
        public static IDictionary<Type, Tuple<object, object>> register;

        public static void Register<Point>(get_x<Point> gx, get_y<Point> gy)
        {
            if (register == null)
            {
                register = new Dictionary<Type, Tuple<object, object>>();
            }
            register[typeof(Point)] = new Tuple<object, object>(gx, gy);
        }
    }


    // The generic library, with functions as "Distance", working on IPoint or on any point
    class Algorithm
    {
        // (1) Using an interface to access coordinates
        public static double Distance(IPoint p1, IPoint p2)
        {
            double dx = p1.x() - p2.x();
            double dy = p1.y() - p2.y();
            return Math.Sqrt(dx * dx + dy * dy);
        }

        // (7) use delegates to access coordinates of not interfaced point type
        public static double Distance<Point>(Point p1, Point p2,
            get_x<Point> gx,
            get_y<Point> gy)
        {
            double dx = gx(p1) - gx(p2);
            double dy = gy(p1) - gy(p2);
            return Math.Sqrt(dx * dx + dy * dy);
        }

        // (8) use traits/accessor to access coordinates of not interfaced point type
        public static double Distance<Point>(Point p1, Point p2, ITraits<Point> traits)
        {
            double dx = traits.get_x(p1) - traits.get_x(p2);
            double dy = traits.get_y(p1) - traits.get_y(p2);
            return Math.Sqrt(dx * dx + dy * dy);
        }

        // (9) use traits/REGISTER accessor to access coordinates of not interfaced point type
        public static double DistanceTR<Point>(Point p1, Point p2)
        {
            ITraits<Point> traits = TraitsRegister.register[typeof(Point)] as ITraits<Point>;
            double dx = traits.get_x(p1) - traits.get_x(p2);
            double dy = traits.get_y(p1) - traits.get_y(p2);
            return Math.Sqrt(dx * dx + dy * dy);
        }

        // (10) use (tupled) lambda/REGISTER accessor to access coordinates of not interfaced point type
        public static double DistanceLR<Point>(Point p1, Point p2)
        {
            var pair = LambdaRegister.register[typeof(Point)] as Tuple<object, object>;
            get_x<Point> gx = pair.Item1 as get_x<Point>;
            get_y<Point> gy = pair.Item2 as get_y<Point>;
            double dx = gx(p1) - gx(p2);
            double dy = gy(p1) - gy(p2);
            return Math.Sqrt(dx * dx + dy * dy);
        }
    }


    // (1) implement IPoint
    class MyPoint1 : IPoint
    {
        private double m_x, m_y;
        public MyPoint1(double x, double y)
        {
            m_x = x;
            m_y = y;
        }
        public double x() { return m_x; }
        public double y() { return m_y; }
    }


    // Non-interfaced class...
    class MyPoint2
    {
        public double x { get; set; }
        public double y { get; set; }
    }

    // (2) Create a wrapper (implementing IPoint)
    class PointWrapper : IPoint
    {
        MyPoint2 m_point;
        public PointWrapper(MyPoint2 p)
        {
            m_point = p;
        }

        public double x() { return m_point.x; }
        public double y() { return m_point.y; }
    }

    // (3) Create a generic wrapper using delegates, implementing IPoint
    class DelPointWrapper<Point> : IPoint
    {
        Point m_point;
        get_x<Point> m_gx;
        get_y<Point> m_gy;

        public DelPointWrapper(Point p, get_x<Point> gx, get_y<Point> gy)
        {
            m_point = p;
            m_gx = gx;
            m_gy = gy;
        }

        public double x() { return m_gx(m_point); }
        public double y() { return m_gy(m_point); }

        public void increase()
        {
            //m_x += 0.000001;
            //m_y += 0.000001;
        }

    }

    // (7) define delegates in a static class
    static class Accessor
    {
        public static get_x<MyPoint2> getX = p => p.x;
        public static get_y<MyPoint2> getY = p => p.y;
    }


    // (8,9) use a sort of Traits class
    class Traits : ITraits<MyPoint2>
    {
        public double get_x(MyPoint2 p)
        {
            return p.x;
        }
        public double get_y(MyPoint2 p)
        {
            return p.y;
        }
    }

    class Program
    {
        const int n = 10000000;

        static void Timed<Point>(String header, Point a, Point b) where Point : IPoint
        {
            double d = 0.0;
            Stopwatch st = new Stopwatch();
            st.Start();
            for (int i = 0; i < n; i++)
            {
                d += Algorithm.Distance(a, b);
            }
            st.Stop();
            System.Console.WriteLine(header + " " + d + " " + st.Elapsed.ToString());
            System.Threading.Thread.Sleep(500);

        }

        // Creating the delegates
        static double point_get_x(MyPoint2 p)
        {
            return p.x;
        }
        static double point_get_y(MyPoint2 p)
        {
            return p.y;
        }

        static void Main(string[] args)
        {

            int option = args.Length > 0 ? Int32.Parse(args[0]) : 10;
            System.Console.Write("{0} ", option);

            // Point without interface
            MyPoint2 p1 = new MyPoint2();
            p1.x = 3;
            p1.y = 5;
            MyPoint2 p2 = new MyPoint2();
            p2.x = 6;
            p2.y = 9;

            switch (option)
            {
                case 1:
                    {
                      // Option 1, point with interface
                      MyPoint1 ip1 = new MyPoint1(3, 5);
                      MyPoint1 ip2 = new MyPoint1(6, 9);
                      Timed("Using interface", ip1, ip2);
                    } break;
                case 2:
                    {
                      // Option 2a:
                      PointWrapper pw1 = new PointWrapper(p1);
                      PointWrapper pw2 = new PointWrapper(p2);
                      Timed("Using wrapper", pw1, pw2);
                    } break;
                case 3:
                    {
                      // Option 2b:
                      DelPointWrapper<MyPoint2> dpw1 = new DelPointWrapper<MyPoint2>(p1, point_get_x, point_get_y);
                      DelPointWrapper<MyPoint2> dpw2 = new DelPointWrapper<MyPoint2>(p2, point_get_x, point_get_y);
                      Timed("Using delegate-wrapper", dpw1, dpw2);
                    } break;
                case 4:
                    {
                      // Option 2c: with inline delegates
                      DelPointWrapper<MyPoint2> idpw1 = new DelPointWrapper<MyPoint2>(p1
                      , delegate(MyPoint2 p) { return p.x; }
                      , delegate(MyPoint2 p) { return p.y; }
                      );
                      DelPointWrapper<MyPoint2> idpw2 = new DelPointWrapper<MyPoint2>(p2
                      , delegate(MyPoint2 p) { return p.x; }
                      , delegate(MyPoint2 p) { return p.y; }
                      );
                      Timed("Using inline delegate-wrapper", idpw1, idpw2);
                    } break;
                case 5:
                    {
                      DelPointWrapper<MyPoint2> ldpw1 = new DelPointWrapper<MyPoint2>(p1, p => p.x, p => p.y);
                      DelPointWrapper<MyPoint2> ldpw2 = new DelPointWrapper<MyPoint2>(p2, p => p.x, p => p.y);
                      Timed("Using lambda-wrapper", ldpw1, ldpw2);
                    } break;
                case 6:
                    {
                      double d = 0.0;
                      Stopwatch st = new Stopwatch();
                      st.Start();
                      for (int i = 0; i < n; i++)
                      {
                      d += Algorithm.Distance(p1, p2, p => p.x, p => p.y);
                      }
                      st.Stop();
                      System.Console.WriteLine("Using delegate-algorithm+lambda " + d + " " + st.Elapsed.ToString());
                    } break;
                case 7:
                    {
                      double d = 0.0;
                      Stopwatch st = new Stopwatch();
                      st.Start();
                      for (int i = 0; i < n; i++)
                      {
                      d += Algorithm.Distance(p1, p2, Accessor.getX, Accessor.getY);
                      }
                      st.Stop();
                      System.Console.WriteLine("Using accessor met static delegates " + d + " " + st.Elapsed.ToString());
                    } break;
                case 8:
                    {
                      Traits traits = new Traits();

                      double d = 0.0;
                      Stopwatch st = new Stopwatch();
                      st.Start();
                      for (int i = 0; i < n; i++)
                      {
                      d += Algorithm.Distance(p1, p2, traits);
                      }
                      st.Stop();
                      System.Console.WriteLine("Using traits " + d + " " + st.Elapsed.ToString());
                    } break;
                case 9:
                    {
                      Traits traits = new Traits();
                      TraitsRegister.Register<MyPoint2>(traits);

                      double d = 0.0;
                      Stopwatch st = new Stopwatch();
                      st.Start();
                      for (int i = 0; i < n; i++)
                      {
                      d += Algorithm.DistanceTR(p1, p2);
                      }
                      st.Stop();
                      System.Console.WriteLine("Using traits/register " + d + " " + st.Elapsed.ToString());
                    } break;
                case 10:
                    {
                      LambdaRegister.Register<MyPoint2>(p => p.x, p => p.y);

                      double d = 0.0;
                      Stopwatch st = new Stopwatch();
                      st.Start();
                      for (int i = 0; i < n; i++)
                      {
                      d += Algorithm.DistanceLR(p1, p2);
                      }
                      st.Stop();
                      System.Console.WriteLine("Using lamba/register " + d + " " + st.Elapsed.ToString());
                    } break;
            }

            System.Console.WriteLine("");
        }
    }
}

And here are the complete timings (for each the smallest of 5 runs):
1 Using interface 50000000 00:00:00.1000114
2 Using wrapper 50000000 00:00:00.1055539
3 Using delegate-wrapper 50000000 00:00:00.1812642
4 Using inline delegate-wrapper 50000000 00:00:00.1843941
5 Using lambda-wrapper 50000000 00:00:00.1810666
6 Using delegate-algorithm+lambda 50000000 00:00:00.1509385
7 Using accessor met static delegates 50000000 00:00:00.1274787
8 Using traits 50000000 00:00:00.1327729
9 Using traits/register 50000000 00:00:00.4192440
10 Using lamba/register 50000000 00:00:00.4194542

Conclusions

Are already done. Specializations in C# are possible. Use them, it makes C# and generics more fun.