含微室

Mimic the interface of golang -- Part 1

I'm learning golang. As a C/C++ programmer, I like the feature of open method and dynamic interface, which are not supported(or not supported in native) by C/C++. In my experience, dynamic interface damage performance a lot, so i write a small bench code:

package main
import (
"fmt"
"time"
)
var (
count = 0
)
func benchmarkFunc(N int, f func(int)) float64 {
start := time.Now()
f(N)
return float64(time.Since(start).Nanoseconds()) / float64(N)
}
func benchmarkFuncPrint(name string, N int, f func(int)) {
ns := benchmarkFunc(N, f)
fmt.Printf("%s %f nano-seconds/call\n", name, ns)
fmt.Printf("%s %f #/second\n", name, 1.e9/ns)
}
type Showable interface {
Show()
}
type Show1 struct {
}
func (p *Show1) Show() {
count++
}
type Show2 struct {
}
func (p *Show2) Show() {
count += 2
}
func call(p interface{}) {
p.(Showable).Show()
}
func call_show(s Showable) {
s.Show()
}
func dynamic(p interface{}, N int) {
for i := 0; i < N; i++ {
// call(p)
p.(Showable).Show()
}
}
func static(s Showable, N int) {
for i := 0; i < N; i++ {
// call_show(s)
s.Show()
}
}
func main() {
N := 10000000
benchmarkFuncPrint("dynamic", N,
func(N int) {
d := Show1{}
dynamic(&d, N)
})
benchmarkFuncPrint("static", N,
func(N int) {
d := Show2{}
static(&d, N)
})
}

This gives

dynamic 19.478761 nano-seconds/call
dynamic 51337966.752537 #/second
static 2.507811 nano-seconds/call
static 398754196.190095 #/second

Which is very quick. Only 19 nano seconds, in my last benchmark, only mutex will cost about 10 nano seconds. And the static version is almost the same as virtual function in C++.

I wrote a toy version of open method for C++, using std::unorderd_map, std::map, google::dense_hash_map, as the storage from std::type_info* to function pointer. And this's the performance counters

lambda: 0.795678 nano-seconds/call
lambda: 1.25679e+09 #/second
unordered_map: 30.4666 nano-seconds/call
unordered_map: 3.28228e+07 #/second
dense hash map: 22.4284 nano-seconds/call
dense hash map: 4.45863e+07 #/second
map: 22.5974 nano-seconds/call
map: 4.42528e+07 #/second

std::map is quicker in this small test, but in larger project, we should divide the performance counter by log(N), where N is the number of function pointers.

It seems that golang's runtime is quicker than google::dense_hash_map + std::mutex. So I felt interesting about it's implementation, and checked the file src/pkg/runtime/iface.c. golang uses a static size(1009) hash table, and only lock if finding first time failed. Chain list in every bucket of hash table is a multiple readers and single writer forward list, reader is wait free.

It's a good implementation, but I don't like the fix size to 1009?(consider a larger project, your performance counter should be divided by ceil(N / 1009), where N is the number of interfaces. And, I guess, using a bit modulo with quadratic open addressing may improve the performance.

So I tried, and I think it did(only test in C++).

template
struct ValueFree {
inline void operator()(const T & v) const {}
};
template <>
struct ValueFree {
inline void operator()(void (*v)(void)) const { }
};

template
struct ValueFree {
inline void operator()(T * v) const {
delete v;
}
};
template
struct TableHash : std::hash {};
template
struct TableHash {
inline size_t operator()(const T*ptr) const {
return (intptr_t)(ptr) >> 4;
}
};
template ,
class DeleteValue=ValueFree>
struct Table {
struct State {
State * old;
size_t size;
size_t buckets;
std::pair table[0];
};
std::atomic state_;
std::mutex add_mutex_;
std::atomic num_find_;
Table() {
const size_t buckets = 8;
auto state = (State*)malloc(sizeof(State) + sizeof(state_.load()->table[0]) * buckets);
state->old = nullptr;
state->size = 0;
state->buckets = buckets;
state_.store(state);
num_find_.store(0);
}
~Table() {
DeleteValue delete_value;
auto state = state_.load();
for (size_t i = 0; i < state->buckets; ++i) {
auto &kv = state->table[i];
if (kv.first != Key()) {
delete_value(kv.second);
}
}
while (state) {
auto old = state->old;
free(state);
state = old;
}
}
Value Find(const Key &k) {
auto state = state_.load(std::memory_order_relaxed);
size_t b = Hash()(k) & (state->buckets - 1);
size_t idx = 0;
while (1) {
auto &o = state->table[b];
if (o.first == k) {
return o.second;
}
if (o.first == Key()) break;
b = (b + ++idx) & (state->buckets - 1);
}
return Value();
}
void Resize(int dir) {
assert(dir == 1);
auto old = state_.load(std::memory_order_acquire);
size_t buckets = old->buckets * 2;
State *state = (State*)malloc(sizeof(State) + sizeof(old->table[0]) * buckets);
state->old = old;
state->size = old->size;
state->buckets = buckets;
for (size_t i = 0; i < buckets; ++i) {
state->table[i].first = Key();
}
for (size_t i = 0; i < old->buckets; ++i) {
auto &o = old->table[i];
if (o.first == Key()) continue;
size_t b = Hash()(o.first) & (state->buckets - 1);
size_t idx = 0;
while (state->table[b].first != Key()) {
b = (b + ++idx) & (state->buckets - 1);
}
state->table[b] = o;
}
state_.store(state, std::memory_order_release);
}
std::pair Add(const Key &k, const Value &value) {
std::lock_guard lk(add_mutex_);
auto state = state_.load(std::memory_order_acquire);
if (state->size * 5 >= 4 * state->buckets) {
Resize(1);
state = state_.load(std::memory_order_relaxed);
}
size_t b = Hash()(k) & (state->buckets - 1);
size_t idx = 0;
while (1) {
auto &o = state->table[b];
if (o.first == k) {
return std::make_pair(o.second, false);
}
if (o.first == Key()) {
o.second = value;
o.first = k;
++state->size;
break;
}
b = (b + ++idx) & (state->buckets - 1);
}
return std::make_pair(value, true);
}
};

Using this as the store to replace google::dense_hash_map, gives:

table: 3.96917 nano-seconds/call
table: 2.51942e+08 #/second

The Table implementation achieve the performance in cost of memory, about half bytes are wasted, and assuming that the value is never changed, deleted. And it also assumes Key() is the empty key.

TODO: open multiple method, with ambiguous resolutions.

reorder to improve cache locality

In large machine learning problems, such as parameter server, we have to calculate over samples. And for every sample, we have to fetch some variables from remote server. To improve performance, we can cache the recent used variables in local memory, and only fetch ones not cached. The problem is that, can we reorder the samples, to minimize the fetching times, and number of fetching variables?

Another case is to reorder a sparse matrix to improve the performance of matrix and vector multiplication, esp. if the sparse matrix is unchanged(or, only values are changed, the sparse structure is kept) and used many times.

The precise problem corresponding to the best cache algorithm. I think is very hard to solve.

A greedy solution is choose the sample from not calculated samples, whose number of fetching variables(not in cache) is smallest.

In the matrix-vector example, if the sparse matrix is not randomly generated, such as user-book rating, there is a lots of people rating almost the same books, then, I guess, the performance improvement may be considerable.

In the parameter server case, the greedy algorithm works for one machine. But we still needs a algorithm to distribute the samples over machines.

代码可读性

不少程序员信奉代码可读性教条, 认为, 代码会被很多人读, 只是偶尔执行.

我的观念恰好相反. 代码会被执行很多遍, 不然, 就不必要写. 被执行很多遍的代码, 需要优化, 极尽能事地优化. 代码量大了时, 要将代码分块:

块间接口要清晰, 块内实现要高效.

被优化的代码, 可能你自己也将看不懂. 那没关系, 你只要能弄明白每一块代码的意图, 而不需要时时让每个人清楚这块代码每行的原理. 不少程序员只是机器, 会解码每行代码, 却不懂合起来的意思, 还责怪别人的代码难懂. 给定代码意图, 只写让所有程序员都理解的代码, 无谓地低效.

如果某人不清楚某块代码的意图, 那是注释的事, 你需要加注释, 以防将来你自己也忘了. 如果某人不清楚某块代码的实现原理, 那是他自己的事. 机器也不懂, 它不需要懂, 它只需要执行. 不懂某块代码的实现原理, 但知道意图, 就很容易重写, 维护问题也就不存在.

自志, 请勿对号入座.

附四年前的代码一段. 估计大部分人不会(立即)清楚这个实现原理, 但清楚意图(仅名字 Single value Linear Regression Analysis 就够了). 当然, 这段还不够高效.

class SLRA:
    '''
    SLRA: Single value Linear Regression Analysis
    usage:
      # Initial a SLRA:
      slra = SLRA(dim)
      # For every incoming vector x of dimension dim, corresponding
      # scalar result y, you can add it to the system by:
      slra.add(x, y)
      # To get the coefficients:
      slra.get()
    '''

    def __init__(self, dim, bigo=1.0e12, eps=1.0e-12):
        self.dim = dim
        self.P = [[i == j and bigo or 0.0 for i in range(dim)] for j in range(dim)]
        self.theta = [eps for i in range(dim)]

    def add(self, x, y):
        px = [sum([self.P[i][j] * x[j] for j in range(self.dim)]) for i in range(self.dim)]
        xpx1 = 1.0 + sum([x[i] * px[i] for i in range(self.dim)])
        K = [px[i] / xpx1 for i in range(self.dim)]
        delta = y - sum([x[i] * self.theta[i] for i in range(self.dim)])
        for i in range(self.dim):
            self.theta[i] += K[i] * delta
            for j in range(self.dim):
                self.P[i][j] -= K[i]*K[j]* xpx1

    def get(self):
        return self.theta

关于 std::tuple 的重排

C++11 的标准里对 tuple 元素的布局没有要求, 但 libstdc++ 的实现里是直接按顺序存放的, 这在部分情况下, 内存布局不太好.

比如, std::tuple<char, double, char> 在 sizeof(double) = 8, alignof(double) = 8 时, 需要占用 24 个字节, 但如果重排为 double, char, char 时, 只需要 16 个字节. 再如 std::tuple<char, int, char, double, char> 需要占用 32 个字节(假设 int 占 4 字节且对齐, double 占 8 字节且对齐), 但重排一下, 可以只占用 16 字节.

重排的 tuple, 需要加一层编译期映射, 保持 std::get 的正确性, 就可以完全替换掉 std::tuple 的实现了.

代码见 gist

blogger compose is too suck to embed the code!

CRTP Pattern with Most Derived Rebind, streambuf as Example

1 The Name

CRTP(Curiously Recurring Template Pattern) is a C++ idiom in which a class X derives from a class template instantiation using X itself as template argument.

CRTP is used for static polymorphism, to avoid virtual functions, and give the compiler a chance to do the best optimization. Call of virtual function doesn't cost much, esp. with correctly predicted branch. but virtual function stop the inline optimization, and cost very much for floating point numeric algorithms.

CRTP is used for equipping object with interfaces based on a very smaller set of interfaces.

2 streambuf Example

std::streambuf is a suitable example. Most of the public/protected functions is implemented based on virtual functions overflow, underflow, uflow. Functions like xsgetn and xsputn with default implementation, based on overflow, underflow, uflow. The pointer gptr, pptr are null if it's not modified by the derived object. Every time in the loop, xsgetn / xsputn check the pointer range. If it's not empty(the pointer range modified by the derived object), memcpy can be used for performance. If it's empty(or all exhausted), uflow / overflow is called.

In some instances(char oriented), gptr and pptr are always null, we can override the implementation without the check. And in some instances, the whole buffer is contiguous, then overflow / uflow always return eof. We can override the virtual functions with these conditions. But not the non-virtual one.

CRTP is suit to solve these two problems, without dynamic cost.

template <class Derived, class Char, class Traits>
struct streambuf_base {
  // ....

  Derived & get_derived() {
    return *this;
  }

  const Derived & get_derived() const {
    return *this;
  }

  int_type snextc() {
    int_type r = traits_type::eof();
    if (get_derived().sbumpc() != traits_type::eof())
      r = get_derived().sgetc();
    return r;
  }

  // ....
};

The default snextc implementation uses sbumpc and sgetc, with get_dervied prefix, the functions are resolved from derived to template base. Only if the derived doesn't override these functions, the template base's default are taken.

For example, stringbuf, alike to std::stringbuf, can implement snextc, sbumpc and sgetc, just with one pointer comparation. All other default implementations will use this new-implemented version, thanks to get_derived prefix.

3 The Problem

The above is fine with just static polymorphism. The problems raise if we want the streambuf_base be compatible with std::streambuf, to reuse all routines who only accepts std::streambuf. So we must inherit from std::streambuf.

We can't inherit from std::streambuf in derived. If we do, any functions not re-implemented by derived will be ambiguous.

Inherit from std::streambuf in streambuf_base is fine. But then the default implementations in streambuf_base will override all in std::streambuf. It's ok in these case, but not acceptable if we inherit from another base, which provides more resonable implementations than streambuf_base.

So, what we want is an order of functions resolution. from derived to ancestor(std::streambuf in former case), and at last streambuf_base.

4 The Solution

First, we must modify the streambuf_base to take another Ancestor type parameters, and provides dummy as default.

We can check whether Ancestor is dummy in every streambuf_base's function. If Ancestor is not dummy, we just call Ancestor's version. Otherwise, the default implementation is taken.

template <class Derived, class Char, class Traits, class Ancestor=dummy>
struct streambuf_base : Ancestor {
  static const bool is_dummy = std::is_same<Ancestor, dummy>::value;
  typedef typename std::conditional<
    is_dummy, streambuf_base, Ancestor>::type ancestor_type;

  int_type snextc() {
    if (!is_dummy) return ancestor_type::snextc();
    int_type r = traits_type::eof();
    if (get_derived().sbumpc() != traits_type::eof())
      r = get_derived().sgetc();
    return r;
  }
};

This introduces some other problems. If Ancestor is another one derived from streambuf_base, the Ancestor::get_derived is wrong. So if Ancestor::snextcis taken, Ancestor::snextc will call another sbumpc and sgetc, not the most-derived object's one.

One solution to this problem is using virtual functions. But make all functions be virutal is not a good solution.

Another solution is to redefine the get_derived. Add another MostDerived paramter to streambuf_base, and add most_derived_rebind meta functor. For every Ancestor, we check whether Ancestor::most_derived_rebind exists. If it exists, streambuf_base inherited from Ancestor::most_derived_rebind::template apply::type instead of Ancestor.

template <class Derived,
          class Ancestor,
          bool=has_nested_type_most_derived_rebind<Ancestor>::value>
struct streambuf_base_rebind_ancestor {
  typedef Base type;
};

template <class Derived, class Base>
struct streambuf_base_rebind_ancestor<Derived, Ancestor, true> {
  typedef typename Ancestor::most_derived_rebind::template apply<Derived>::type type;
};

template <class Derived, class Char, class Traits,
  class Ancestor=dummy,
  class MostDerived=Derived>
struct streambuf_base : streambuf_base_rebind_ancestor<MostDerived, Ancestor>::type {
  typedef typename streambuf_base_rebind_ancestor<
    MostDerived, Ancestor>::type _ancestor_type;
  typedef typename std::conditional<
    std::is_same<_ancestor_type, dummy>::value,
    streambuf_base, _ancestor_type>::type ancestor_type;

  MostDerived & get_derived() {
    return *this;
  }

  const MostDerived & get_derived() const {
    return *this;
  }

  struct most_derived_rebind {
    template <class T>
    struct apply {
      typedef streambuf_base<Derived, Char, Traits, Ancestor, T> type;
    };
  };
};

5 The Result

With above streambuf_base, we can implements text format input/output routines in very higher performance, even faster than strtod / strtol in experiments(locale is ignored).

The floating pointer converion is done by double-conversion, in shorted mode. For printf/scanf and stringstream, the precision is 20.

The random sequence for floating is generated by random() / M_PI and for long, just random().

For fscanf and fprintf, the FILE* is obtained by fmemopen and open_memstream.

Random generated input speed(in nano-seconds)
implements / type	long	float	double
streambuf_base	45.13806234741213	87.06188110351565	193.94798608398446
strtox	56.38690736389162	99.91599140930177	246.5969282989502
fscanf	201.75892387390138	244.81490730285645	437.06397615051276
istringstream	168.8839549407959	382.54306460571297	674.1460767364503

Random generated output speed(in nano-seconds)
implements / type	long	float	double
streambuf_base	109.96006123352053	221.0190715637207	285.9951219940186
snprintf	290.8328487548828	1263.6021590270996	1915.7261646881104
fprintf	279.30502569580074	1298.2731170043944	2073.194082229614
ostringstream	122.75186352539063	1086.216874282837	1651.6398792114258

simple/naive implementation of continuation/coroutine/generator on GNU/Linux

1 continuation

In some senses, ucontext.h is already a implementation of continuation, which is available on most modern OS(standarded by POSIX). In C99, we also have setjmp/longjmp to implement continuation. But setjmp/longjmp is weaker than ucontext.h.
Folowing code are just a demo, not suit for production usage. g++-4.7 is required.

#include 
#include  // for type erasure, or, we have to make continuation a template

#define CLOG(...) // we don't have/need CLOG

struct continuation {
  template <class Func>
  continuation(const Func & func, size_t stack_size=1<<18) {
    m_stack.resize(stack_size);
    m_func = func;
    int r = getcontext(&m_ctxs[1]);
    assert(r == 0);
    m_ctxs[1].uc_link = &m_ctxs[0];
    m_ctxs[1].uc_stack.ss_sp = &m_stack[0];
    m_ctxs[1].uc_stack.ss_size = m_stack.size();
    // TODO: x86_64 is safe to pass pointer with glibc, but for other platform ...
    makecontext(&m_ctxs[1], (void (*)(void))&s_context_func, 1, this);
  }

  continuation(const continuation &r) = delete;
  continuation & operator=(const continuation &r) = delete;

  int swap() {
    if (m_in_done) return -1;
    CLOG(DEBUG, "swap from", (m_is_in ? "in" : "out"));
    m_is_in = !m_is_in;
    swapcontext(&m_ctxs[!m_is_in], &m_ctxs[m_is_in]);
    CLOG(DEBUG, "after swap");
    if (m_in_done) return 0;
    return 1;
  }

  static void s_context_func(continuation *self) {
    return self->context_func();
  }

  void context_func() {
    m_func();
    m_in_done = 1;
  }

  ucontext_t m_ctxs[2];
  std::function<void(void)> m_func;
  std::vector<char> m_stack;
  bool m_is_in = false;
  int m_in_done = 0;
};

2 generator

template <class T>
struct generator : continuation {
  template <class Func>
  generator(const Func &func)
  : continuation(func) {

  }

  void yield(const T &value) {
    m_value = value;
    swap();
  }

  T yield() {
    if (swap() <= 0) throw std::runtime_error("end of iteration");
    return m_value;
  }

  T volatile m_value;
};

3 coroutine

struct coroutine : continuation {
  template <class Func>
  generator(const Func &func)
  : continuation(func) {

  }

  void yield(coroutine &other) {
    other.swap();
  }
};

4 usage and test

4.1 primes

A simple program to print primes(not all).

bool is_prime(int n) {
  if (n % 2 == 0) return n == 2;
  for (int p = 3; p * p <= n; p += 2) {
    if (n % p == 0) return false;
  }
  return true;
}

void primes(generator<int> &g) {
  g.yield(2);
  int i = 3;
  while (1) {
    if (is_prime(i)) g.yield(i);
    i += 2;
  }
}

int main() {
  generator<int> g([&]() {
      primes(g);
    });
  while (1) {
    int p = g.yield();
    std::cout << p << std::endl;
  }
}

4.2 producer and consumer

void producer(int n, int *x, continuation & cont) {
  for (int i = 0; i < n; ++i) {
    *x = i * i;
    cont.swap();
  }
}

void consumer(int volatile *x, continuation &cont) {
  while (cont.swap() > 0) {
    CLOG(INFO, "got", *x);
  }
}

向量和矩阵的表达式模板

不是我折腾矩阵和向量, 是矩阵和向量折腾我.
早期为速度, 使用 cblas 比较多. cblas 非常不好用, double, float 函数各一还好说, 关键是大部分函数的参数太多了, 如

void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
                 const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
                 const int K, const double alpha, const double *A,
                 const int lda, const double *B, const int ldb,
                 const double beta, double *C, const int ldc);

基本上, 每次使用, 我都测试一下 M, N, K 几个值差得比较大的情形, 看看有没有段错误, 这个方法比认真检查参数还要管用.
后来, 用 C++ 比较多, 仿照 STL, 希望尽量泛化, 不要 double, float 各写一遍. 还希望进一步简化, 对所有矩阵, dense 的也好, csr, csc 也好, 我们应该只需要写 gemm(A, B).
还要将计算和存储分开, 比如 gemm(A, B) 只表示一个表达式, 并不做任何实际运算, 需要时, 可以将结果并行地存到另一个区域. 这样, 比如写 C=gemm(A, B) * alpha + beta * C 就比上面的 cblas 简明易懂多了. 还可以用 C.parallel_assign(gemm(A, B) * alpha + beta * C) 表示并行地计算
$$C = \alpha A B + \beta C.$$
要想实现上面的表达式, 我们需要一个 proxy 抽象, vector_proxy, matrix_proxy. 矩阵, 向量的所有公有操作, 在 proxy 里实现. 具体的向量, 矩阵继承此 proxy 即可.
对所有 vector_proxy, 定义

template <class Derived, class Value, class Reference>
struct vector_proxy {
...
size_t size() const;
value_type get(size_t idx) const;
reference_type get(size_t idx);
...
};

其中, Derived 为具体的向量, 需要定义 get, size, 并继承 vector_proxy<Derived>.
再为 vector_proxy 定义加, 减等基本运算, 结果也为一个 proxy, 此 proxy 的 get(size_t idx) 函数实现计算.
从 proxy 到 proxy 的赋值如下.

template <class Derived, class Value, class Reference>
struct vector_proxy {
...
 
template <class R>
void assign(const vector_proxy<R> &r) {
  size_t n = size();
  assert(r.size() == n);
  for (size_t i = 0; i < n; ++i) {
     get(i) = r.get(i);
  }
}
 
template <class R>
void parallel_assign(const vector_proxy<R> &r) {
  size_t n = size();
  assert(r.size() == n);
#pragma omp parallel for 
  for (size_t i = 0; i < n; ++i) {
    get(i) = r.get(i);
  }
}
 
...
 
}

matrix_proxy 和 vector_proxy 完全类似.
具体到 gemm_proxy 上, 如下实现 get 方法即可.

value_type get(size_t i, size_t j) {
  return vec_inner_prod(m_A.row(i), m_B.col(j));
}

矩阵的 row, col 返回的也是 vector_proxy. 在 vec_inner_prod 里, 根据向量是 dense 还是 sparse, 采用不同的计算方法即可.
gemm_proxy 里的 row(i), col(j) 可以采用矩阵和向量的乘积 gemv(m_B, m_A.row(i)), gemv(m_A, m_B.col(j)). 在 gemv_proxy 里, 还是用 vec_inner_prod 来自动区分 dense/sparse 的向量.
gemm(A, B) 在 A, B 都是稀疏矩阵时, 情况比较复杂, 使 gemv 返回一个稀疏向量的实现效率比较低, 我直接将结果当成 dense 了.