Bannalia: trivial notes on themes diverse

60/40 split

2009-06-15T22:20:00.002+02:00

Tim Harford recounts in his book The Logic of Life an anecdote about the authors of Freakonomics, Steven D. Levitt and Stephen J. Dubner: when discussing the terms of their collaboration for writing the book, Levitt stated that he wouldn't go for less than a 60/40 split on revenues; Dubner, on the other hand, would settle for no less than 60/40 either. The apparent confrontation vanished when they realized that each was meaning the other to get 60%, so they decided to join the venture.

Charming as it may sound, this story cannot be true: if I'm settling for 40% of the revenues and my partner proposes a 60/40 split, the proposal immediately fits my goals no matter whom I assume will get each part of the split.

About break statements

2009-03-23T04:49:00.000+01:00

In a recent article, Andrew Koenig argues that break statements are harmful because they make it harder to deduce the program context at loop termination. For instance, given:

while (n != 0) { /* do something */ }

n can be assumed to be zero at loop termination unless there is a break within it. Koenig proposes two different approaches to alleviating this problem:

Include the program context at the break statement as a possible value for the program context at loop termination.
Force the loop termination condition just before breaking.

There is another alternative: exclude the loop termination condition from the evaluable program context. Here is a way to do that:

for (int local_n = n; local_n != 0;) { /* do something with local_n */ }

As the scope of local_n does not extend outside the loop, program context after loop termination does not depend on the termination procedure. A cruder, but also more flexible way to do the same is the following:

‎{
‎  int local_n = n;
‎  while (local_n != 0) { /* do something with local_n */ }
‎}

The somewhat unexpected surrounding braces are a conspicuous reminder that local_n does not form part of the global program state.

A curious syntactic transformation

2009-03-17T18:01:00.000+01:00

Consider the sentence

Rationalism has many followers.

What in principle looks like a rather dull N+V+O statement does not however allow for seeminlgy innocuous variations on the object:

*Rationalism has many people.

The reason why the former is valid while the latter is not is that followers are followers of something, in this particular case rationalism:

Rationalism has many followers [of rationalism].

In fact, we can view our sentence as a mere rewording of:

There are many followers of rationalism,

which leads us to hypothetize the existence of the following syntactic transformation:

There is/are [Det] N of NP → NP has/have [Det] N.

For instance:

There were lots of fans of the Beatles → The Beatles had lots of fans.
There are no enemies of Rome → Rome has no enemies.

This transformational rule explains from a purely syntactical perspective why superficially similar sentences like *Rationalism has many people are invalid --they have no "There is..." equivalent.

Interestingly, the rule seems to operate in other languages apart from English (I presume that at least in most modern European languages):

El racionalismo tiene muchos seguidores.
Le rationalisme a beaucoup d'adeptes.
Rationalismus hat viele Anhänger.

This points to some general (maybe universal?) mechanism of reification of the possesion relationship between nouns.

The mythical Bell Curve in Human Resources

2009-03-10T19:50:00.002+01:00

Suppose a company's HR department wishes to establish a statistical model for their personnel performance in order to set up outcome predictions for the company's bonus system. It is all too easy to assume that performance will follow a normal distribution such as this:

Fig. 1: Performance distribution under the Bell Curve assumption.

The assumption stems from the deeply held custom in Psychology and Sociology of using the Bell Curve to model any a priori unknown human trait. Actually, if the company's hiring process is indeed efficient in selecting better than average people, this assumption is a complete contradiction.

Consider a simplistic hiring process in which performance is assesed by means of a test, so that applicants scoring some fixed minimum at the test are hired, and let us concede for the sake of the argument that the test is a perfect predictor of performance. The resulting performance distribution is depicted at the figure.

Fig. 2: Performance distribution with a test-based selection process.

The distribution has positive skewness, i.e. its right tail is longer than the left tail. So, the normal approximation with the same mean and variance (shown in dotted line) both overestimates low performers and underestimates high performers. It also underestimates low-to-normal performers and overestimates normal-to-high performers.

In other hiring process scenario, the best among N candidates is selected. The resulting distribution is depicted at the following figure for N = 10.

Fig. 3: Performance distribution with best-of-10 selection process.

We have positive skewness again, though not as marked as in the prvious case. Skewness grows as N does. Again, the normal approximation results in overestimation of low performers and underestimation of high performers.

Finally, we consider a two-stage hiring process where applicants are first filtered by a test and then the best candidate out of N is selected.

Fig. 4: Performance distribution with test prefiltering and best-of-10 selection process.

The test filter results in a slightly larger positive skewness. As in previous cases, normal approximation predicts more low performers and less high performers than the real case.

To summarize: hiring process not only results in a personnel performance distribution with a higher than average mean (which is the primary purpose of any hiring process); the distribution will also have positive skewness, with more excellent and less deficient people than predicted by the Bell Curve.

More on price gas hysteresis

2008-11-28T05:47:00.006+01:00

We redo the hysteresis analysis we conducted on gas prices, now for a longer period of time and using much more data than before:

Average weekly retail prices without taxes of 95 octane gasoline and gasoil in Spain have been obtained from the European Commision Oil Bulletin (thanks to Wonka for the pointer).
The US Energy Information Administration provides daily spot prices of Brent oil in dollars per barrel.
Historical data on euro to dollar exchange rates can be retrieved from FXHistory.

Using this information we have constructed the following figure showing the evolution of Brent prices and gasoline and gasoil Spanish retail prices without taxes, all in c€ per liter, between January 2006 and October 2008.

And this is the scatter plot of Δ(gasoline price before taxes) against Δ(Brent price):

The linear regressions of this plot for the semiplanes Δ(Brent price) ≥ 0 and ≤ 0 are:

ΔBrent ≥ 0 → y = f₊(x) = b₊ + m₊x = 0.1854 + 0.1691x,
ΔBrent ≤ 0 → y = f₋(x) = b₋ + m₋x = 0.1273 + 0.3636x.

On one hand, m₋ is more than twice as large as m₊, which indicates that gasoline price follows oil price reductions much more strongly than oil price increases when these variations are large; when oil price variations are small, it is the fact that both b₊ and b₊ are positive that dominates the situation: gasoline prices have a very clear bias towards increasing then. The "consumer unfavorable" region in which f₊(x) > −f₋(−x) is

−1.61 c€/l < ΔBrent < 1.61 c€/l.

As for gasoil, the plot look like this:

with regression lines

ΔBrent ≥ 0 → y = f₊(x) = b₊ + m₊x = 0.2669 + 0.1269x,
ΔBrent ≤ 0 → y = f₋(x) = b₋ + m₋x = 0.2391 + 0.4134x.

Again, m₋ is much greater than m₊ (~3.25 times), but this is masked for small variations of Brent price by the constants b₊ and b₊ being greater than zero. The consumer unfavorable region is

−1.77 c€/l < ΔBrent < 1.77 c€/l,

which is very similar to the consumer unfavorable region for gasoline.

So, we have a somewhat mixed scenario: when oil prices do not vary much, there is a tendency for retail gas prices to increase, whereas when oil price variations are high, reductions are tracked much more strongly than increases. We propose two explanations, both highly speculative, for this phenomenon:

Small variations of oil prices do not reach public attention, which allows retailers to silently increase gas prices; high variations in oil price, however, generally make it to the news, so that retailers are under greater public scrutiny when transferring these variations to the pump.
So as to soften the public negative reaction to price increments (and associated contractions in consumption), retailers do not reflect high increases in oil price as abruptly as oil price reductions, but opt to distribute the increment over a longer period of time; this would account both for the fact that m₊ < m₋ and for the positive bias in small price variations.

Optimizing red-black tree color bits

2008-11-24T21:31:00.000+01:00

C++ associative containers are almost universally implemented on top of a structure known as a red-black tree. An rb tree node has approximately the following look:

‎enum color_type{red=false,black=true};
‎
‎template<typename T>
‎struct rb_node
‎{
‎  color_type color;
‎  rb_node*   parent;
‎  rb_node*   left;
‎  rb_node*   right;
‎  T          value;
‎};

The color of a node, from which the data structure receives its name, is used to implement some algorithms to keep the tree balanced. Note that the information conveyed by color is exactly one bit; yet, this member takes at least one char of storage and, in most platforms, much more than that because of alignment issues. Increasingly more libraries implement an optimization technique that is able to get rid of the extra storage demanded by color by embedding the information bit into one of the additional node fields. As this technique does not appear to be much documented on the Internet, this entry describes it in detail. Much of the material draws from the source code of Boost.MultiIndex. Some acquantaince with Boost.MPL is welcome to understand what follows.

To embed a bit into the representation of a pointer we need to use in place of the pointer an unsigned integer type with the same size. In most 32bit platforms this type is unsigned int, but the standard does not guarantee this, or even that such an integral type exists. We assume that we have the following utilities at our disposal:

‎typedef ... has_uintptr_type;
‎typedef ... uintptr_type;

where has_uintptr_type is a Boost.MPL boolean indicating whether an unsigned integral type exists with exactly the same size as a pointer, and if so uintptr_type is such type. We will provide later an implementation for these.

Alignment considerations dictate that many types must lie in memory addresses that are multiple of some small integer greater than 1; in particular, it is extremely usual that the alignment of a type T is even (i.e. a multiple of 2), in which case the least significant bit of the representation of pointers to T is always zero. And this is precisely the situation in which we can reuse this bit to embed the color information:

‎template<typename T,bool ColorBitCompression=true>
‎class rb_node:
‎  public boost::mpl::if_c<
‎    ColorBitCompression&&
‎    has_uintptr_type::value&&
‎    (boost::alignment_of<compressed_rb_node_header>::value%2==0),
‎    compressed_rb_node_header,
‎    regular_rb_node_header
‎  >::type
‎{
‎  T value_;
‎
‎public:
‎  T&       value(){return value_;}
‎  const T& value()const{return value_;}
‎};

rb_node uses a header with the color bit impression technique when:

This is requested by the user via ColorBitCompression.
There is an unsigned integer with which pointers can be represented.
The alignment of headers is even.

Otherwise a regular header is used. The implementation of the latter is entirely obvious:

‎class regular_rb_node_header
‎{
‎  typedef regular_rb_node_header* pointer;
‎
‎  color_type                color_;
‎  regular_rb_node_header*   parent_;
‎  regular_rb_node_header*   left_;
‎  regular_rb_node_header*   right_;
‎
‎public:
‎  color_type& color(){return color_;}
‎  color_type  color()const{return color_;}
‎  pointer&    parent(){return parent_;}
‎  pointer     parent()const{return parent_;}
‎  pointer&    left(){return left_;}
‎  pointer     left()const{return left_;}
‎  pointer&    right(){return right_;}
‎  pointer     right()const{return right_;}
‎};

As for compressed_rb_node_header, the most interesting part lie in the use of proxy classes color_ref and parent_ref to isolate the user from the representation of the color information and the parent pointer as a merged uintptr_type:

‎class compressed_rb_node_header
‎{
‎  typedef compressed_rb_node_header* pointer;
‎  struct color_ref
‎  {
‎    color_ref(uintptr_type* r_):r(r_){}
‎ 
‎    operator color_type()const
‎    {
‎      return color_type(*r&uintptr_type(1));
‎    }
‎ 
‎    color_ref& operator=(color_type c)
‎    {
‎      *r&=~uintptr_type(1);
‎      *r|=uintptr_type(c);
‎      return *this;
‎    }
‎ 
‎    color_ref& operator=(const color_ref& x)
‎    {
‎      return operator=(x.operator color_type());
‎    }
‎ 
‎  private:
‎    uintptr_type* r;
‎  };
‎  struct parent_ref
‎  {
‎    parent_ref(uintptr_type* r_):r(r_){}
‎ 
‎    operator pointer()const
‎    {
‎      return (pointer)(void*)(*r&~uintptr_type(1));
‎    }
‎ 
‎    parent_ref& operator=(pointer p)
‎    {
‎      *r=((uintptr_type)(void*)p)|(*r&uintptr_type(1));
‎      return *this;
‎    }
‎ 
‎    parent_ref& operator=(const parent_ref& x)
‎    {
‎      return operator=(x.operator pointer());
‎    }
‎
‎    pointer operator->()const
‎    {
‎      return operator pointer();
‎    }
‎
‎  private:
‎    uintptr_type* r;
‎  };
‎
‎  uintptr_type parentcolor_;
‎  pointer      left_;
‎  pointer      right_;
‎
‎public:
‎  color_ref  color(){return color_ref(&parentcolor_);}
‎  color_type color()const{return color_type(parentcolor_&std::size_t(1ul));}
‎  parent_ref parent(){return parent_ref(&parentcolor_);}
‎  pointer    parent()const
‎    {return (pointer)(void*)(parentcolor_&~uintptr_type(1));}
‎  pointer& left(){return left_;}
‎  pointer  left()const{return left_;}
‎  pointer& right(){return right_;}
‎  pointer  right()const{return right_;}
‎};

A complete example program is provided. In typical 32bit architectures, color bit compression makes the size of rb_node<T> reduce by 4 bytes.

This color bit compression technique introduces some penalty when accessing and modifying the color and parent information. This penalty, however, is entirely negligible. What is more, the reduction of size brought about by color bit compression usually results in a performance speedup for associative containers taking advantage of this trick: as nodes are smaller, more of them can be stored in L1 cache, which accelerates execution noticeably.

One final remark about this technique: as the parent pointer is no longer declared as a regular pointer, but rather represented by a masked integer, most C++ garbage collectors will not recognize the implicit pointer, which in principle can cause problems. All is fine, however, as it is only the parent that is masked: the right and left pointers are stored as genuine pointers, and as every node of a red/black tree is reachable only through left and right garbage collectors will not fail.

Thank you to Thorsten Ottosen for discussing this color bit compression technique with me.

Infinitely many voters

2008-11-20T18:51:00.000+01:00

As we have discussed before, the famous Arrow's Theorem is equivalent to this basic theorem on ultrafilters:

Every ultrafilter on a finite set is principal.

(See for instance Joel W. Robbin's paper for a proof of this equivalence.) Basically, the set on which the ultrafilter is built represents the voters, and the member sets of an ultrafilter can be equated with forcing coalitions, groups of voters which, if unanimous on a given preference, determine the overall result with respect to that preference. The theorem above says that every ultrafilter on a finite set of voters has a forcing coalition of size 1, i.e. a dictator.

If we have an inifite set, however, the Axiom of Choice implies that there are non-principal ultrafilters: so if there were infinitely many voters we could devise (alas, nonconstructively) a fair voting system free of dictatorships. It is interesting to try to take advantage of this result to somehow circumvent the limitations for the finite case, and see what fails (something must fail after all):

Let our voters be represented by the set V={v₁,...,v_N}: we "extend" them by assigning to each v_i an enumerable set of fictitious voters v_i1, v_i2,... Thus the set V' of fictitious voters is infinite and we can have a non-principal ultrafilter F on it. We design the following voting system: for any given preference the vote of each v_i is obtained and assigned to the associated fictitious voters v_ij; if the overall resulting set of fictitious voters supporting the preference is a forcing coalition according to F, the preference is approved, otherwise rejected.

We already know that this voting system cannot be fair, and in fact it is easy to uncover the underlying flaw: the range of our mapping f : V → V' is not the entire powerset of V', but only a finite subset of it (with cardinality 2^N); although F is non-principal, it becomes principal (streching somewhat the meaning of the concept) when restricted to range(f): F necessarily contains then a pseudodictatorial forcing coalition comprising exactly all the fictitious voters associated to a single v_i, who is the resulting dictator for the so devised voting system.

State semantics

2008-11-16T20:32:00.000+01:00

The following is a proposal for the classification of types in an imperative language according to the way the state of their objects changes over time.

Immutable. Objects of an immutable type cannot change state during their lifetime.

Stateful. On the other side of the spectrum, an object of a stateful type publish some methods that cause the object's internal state to change when executed.

Bound. (For lack of a better name.) The state of an object x of a bound type remains immutable unless x is explicitly given the same state as some y through explicit assignment x = y. A very popular example of a bound type is java.lang.String. Although Java strings are commonly characterized as immutable, it is actually their associated states that are immutable, since an object of type java.lang.String can be rebound through assignment. In C++, bound semantics are approximated by pointers to constant objects—C++ references, on the other hand, do not have bound semantics since they cannot be reassigned, that is, rebound.

It turns out that bound types are sufficiently powerful to emulate stateful semantics using the following mapping technique: let T be a stateful type and an associated mutable method

void T::f(args);

We construct an associated bound type B_T whose internal state is represented by values of T. We replace T::f with the following immutable method:

‎B_T B_T::f(args)
‎{
‎  T state=this->state;
‎  state.f(args);
‎  return B_T(state);
‎}

so that the expression

‎T x;
‎...
‎x.f(args);

can be rewritten with B_T as

‎B_T x;
‎...
‎x=x.f(args);

And similarly we can provide a stateful façade to a bound type: let T be a bound type with some method

‎T T::f(args);

We define a stateful type S_T whose internal state is represented by values of T, and implement the associated method f as follows:

‎void S_T::f(args)
‎{
‎  this->state=this->state.f(args);
‎}

In summary, bound and stateful types are to some extent equivalent, though practical considerations might favor one kind over the other depending on the use context. For instance, bound types are more amenable to interning techniques and can be more easily equipped with locking mechanisms for use in multithreaded environments.

Gas price hysteresis

2008-11-12T08:04:00.001+01:00

Due to the extreme rise of gas prices during the last couple of years, car users have become very sensitive (if not specially elastic) to the translation of oil price variations to retail gasoline prices. In particular, some have expressed the suspicion that increases in oil price are rapidly transferred to gas stations, whereas when the oil price decreases, this reflects less directly into a reduction of retail gas prices. Is this a case of actual price hysteresis or just consumer hysteria?

The figure shows the evolution of Brent prices and average retail prices of 95 octane gasoline and gasoil (diesel) in Spain, all in c€/l (source: Cores, via WonkaPistas).

The final price of a liter of gasoline or gasoil in Spain can be expressed as

final price = (1 + VAT) × (price before taxes + special tax),

where the applicable VAT is currently 16% and the special tax, known as the impuesto especial sobre hidrocarburos, is approximately fixed. So, the variations on this price do not depend on this special tax:

Δ(final price) = (1 + VAT) × Δ(price before taxes).

These variations can be then attributed to the volatile part of the product cost, the most important of which is oil price. This is how Δ(gasoline price before taxes) plots against Δ(Brent price) monthly between February 2007 and August 2008:

The dotted lines are the linear regressions of the points for Δ(Brent price) ≥ 0 and ≤ 0, respectively:

ΔBrent ≥ 0 → y = f₊(x) = b₊ + m₊x = −0.5378 + 0.9885x,
ΔBrent ≤ 0 → y = f₋(x) = b₋ + m₋x = 2.4610 + 1.4865x,

How can we interpret this? If the variations of gasoline price were symmetrical with respect to the sign of Δ(Brent price), we would have

f₊(x) = −f₋(−x),
b₊ = −b₋,
m₊ = m₋,

which is clearly not the case. In fact, we have that

f₊(x) > −f₋(−x) for x in [0,3.86),

that is, variations in the price of the gasoline when 0 ≤ ΔBrent < 3.86 are larger than the corresponding reductions when − 3.86 < ΔBrent ≤ 0. In this interval, gasoline price can be said to present hysteresis unfavorable to the consumer. The situation is reversed when ΔBrent > 3.86, though this phenomenom is mainly accounted for by the presence of the outlier point at (−5.64,−5.41).

As for the variations on the price of gasoil before taxes, we have the following:

with regression lines

ΔBrent ≥ 0 → y = f₊(x) = b₊ + m₊x = −0.7877 + 1.3841x,
ΔBrent ≤ 0 → y = f₋(x) = b₋ + m₋x = 2.4676 + 1.6557x,

resunting in an user unfavorable hysteresis interval −6.19 < ΔBrent < 6.19, which spans the entire dataset.

In conclusion, only to some extent does the analysis corroborate the feeling many users share that oil companies do not reflect oil price reductions as faithfully as they do reflect increases. These results are not very trustworthy since negative variations of oil price are underrepresented. Given the recent drop in the international oil market, it will be interesting to repeat the calculations when the data for the last two months become available.

Multiattribute querying with Boost.MultiIndex

2008-11-08T04:55:00.002+01:00

We saw in a prior entry how to define database indexes on a set of columns so that all possible queries involving clauses of the form column_i = value_i are accelerated by some index. We see now how to model this situation for an in-memory C++ container, which gives us the opportunity to use Boost.MultiIndex and show some template and preprocessor metaprogramming in action.

So, suppose we have the following struct:

‎struct element
‎{
‎  int x1,x2,x3,x4;
‎};

and we are assigned the task of designing a container of elementss along with efficient querying functions of the form:

‎typedef ... container;
‎
‎template<int N>
‎struct field
‎{
‎  field(int x):x(x){}
‎  int x;
‎};
‎
‎template<int N0>
‎const element* find(const container& c,const field<N0>& f0);
‎
‎...
‎
‎template<int N0,int N1,int N2,int N3>
‎const element* find(
‎  const container& c,
‎  const field<N0>& f0,const field<N1>& f1,
‎  const field<N2>& f2,const field<N3>& f2);

so that operations like the following are possible:

‎container c;
‎...
‎// x3==10 && x1==20
‎find(c,field<3>(10),field<1>(20));
‎
‎// x1==5 && x4==40 && x2==11
‎find(c,field<1>(5),field<4>(40),field<2>(11));

Boost.MultiIndex is a natural choice to implement the container. As we already know, we need as many as six indices to cover all the possible variations:

‎typedef member<element,int,&element::x1> key1;
‎typedef member<element,int,&element::x2> key2;
‎typedef member<element,int,&element::x3> key3;
‎typedef member<element,int,&element::x4> key4;
‎
‎‎typedef multi_index_container<
‎  element,
‎  indexed_by<
‎    ordered_non_unique<composite_key<element,key1,key3,key2,key4> >,
‎    ordered_non_unique<composite_key<element,key4,key1,key3,key2> >,
‎    ordered_non_unique<composite_key<element,key2,key1,key3,key4> >,
‎    ordered_non_unique<composite_key<element,key4,key2,key1,key3> >,
‎    ordered_non_unique<composite_key<element,key3,key2,key1,key4> >,
‎    ordered_non_unique<composite_key<element,key4,key3,key2,key1> >
‎  >
‎> container;

Actually, the combined composite keys are redundant, and we can shorten some of them. Also, we add a little instrumentation to the type for future use, so the final definition of container is:

‎struct key1:member<element,int,&element::x1>,boost::mpl::int_<1>{};
‎struct key2:member<element,int,&element::x2>,boost::mpl::int_<2>{};
‎struct key3:member<element,int,&element::x3>,boost::mpl::int_<3>{};
‎struct key4:member<element,int,&element::x4>,boost::mpl::int_<4>{};
‎
‎template<int>struct cover;
‎
‎typedef multi_index_container<
‎  element,
‎  indexed_by<
‎    ordered_non_unique<
‎      tag<cover<1>,cover<13>,cover<123>,cover<1234> >,
‎      composite_key<element,key1,key3,key2,key4>
‎    >,
‎    ordered_non_unique<
‎      tag<cover<4>,cover<14>,cover<134> >,
‎      composite_key<element,key4,key1,key3>
‎    >,
‎    ordered_non_unique<
‎      tag<cover<2>,cover<12> >,
‎      composite_key<element,key2,key1>
‎    >,
‎    ordered_non_unique<
‎      tag<cover<24>,cover<124> >,
‎      composite_key<element,key4,key2,key1>
‎    >,
‎    ordered_non_unique<
‎      tag<cover<3>,cover<23> >,
‎      composite_key<element,key3,key2>
‎    >,
‎    ordered_non_unique<
‎      tag<cover<34>,cover<234> >,
‎      composite_key<element,key4,key3,key2>
‎    >
‎  >
‎> container;

The keyn extractor derives from boost::mpl::int_<n>, which makes these types easily indexable. As for the cover tags, they simply indicate the different query formulas supported by each index using a straightforward coding technique: for instance, the second index covers the querying formulas

‎x4==v4 --> cover<4>
‎x1==v1 && x4==v4 --> cover<14>
‎x1==v1 && x3==v3 && x4==v4 --> cover<134>

To avoid ambiguities, cover numbers are coded assuming the attribute names appear in the order x1, x2, x3, x4, even if the tagged index does not respect this order in the specification of the composite key; also, when two different types cover the same formula, we assign the tag to only one of them. So, there are a total of 15 different tags, exactly as many as the number of different query formulas with up to 4 attributes. cover_number does the metaprogrammatic work of calculating a cover number from an MPL sequence of field<n> types:

‎template<typename FieldSequence>
‎struct cover_number:
‎  boost::mpl::fold<
‎    typename boost::mpl::sort<FieldSequence>::type,
‎    boost::mpl::int_<0>,
‎    boost::mpl::plus<
‎      boost::mpl::times<boost::mpl::_1,boost::mpl::int_<10> >,
‎      boost::mpl::_2
‎    >
‎  >::type
‎{};

The following is a possible realization of the overload of find with three fields:

‎template<int N0,int N1,int N2>
‎const element* find(
‎  const container& c,
‎  const field<N0>& f0,const field<N1>& f1,const field<N2>& f2)
‎{
‎  // find the tag associated to the index conering our fields
‎  typedef cover<
‎    cover_number<
‎      boost::mpl::vector_c<int,N0,N1,N2>
‎    >::value
‎  >                                                tag;
‎
‎  // retrieve the index type
‎  typedef typename container::index<tag>::type     index_type;
‎
‎  // retrieve the type of the associated composite key extractor
‎  typedef typename index_type::key_from_value      composite_key_type;
‎
‎  // an MPL sequence containing the field types passed
‎  typedef boost::mpl::vector<
‎    field<N0>,field<N1>,field<N2>
‎  >                                                fields_type;
‎
‎  // get the index
‎  const index_type&         i=c.get<tag>();
‎
‎  // pack the field values into a tuple...
‎  boost::tuple<int,int,int> t=boost::make_tuple(f0.x,f1.x,f2.x);
‎
‎  // ...and use it with i.find(), rearranging the tuple
‎  // so that each field<n> matches its key<n> in the composite key
‎  typename index_type::iterator it=i.find(
‎    boost::make_tuple(
‎      get<composite_key_type,fields_type,boost::mpl::int_<0> >(t),
‎      get<composite_key_type,fields_type,boost::mpl::int_<1> >(t),
‎      get<composite_key_type,fields_type,boost::mpl::int_<2> >(t)
‎    )
‎  );
‎  if(it!=i.end())return &*it;
‎  else return 0;
‎}

The only missing part is the get template function, that accepts a composite key extractor type, an MPL sequence of fields, a position and a tuple of field values and returns the value associated to the field<n> with the same n as the key<n> at the indicated position in the composite key extractor:

‎template<
‎  typename CompositeKey,typename FieldSequence,
‎  typename Pos,typename Tuple
‎>
‎int get(const Tuple& t)
‎{
‎  typedef typename boost::tuples::element<
‎    Pos::value,
‎    typename CompositeKey::key_extractor_tuple
‎  >::type key_at_pos;
‎  const int M=
‎    boost::mpl::distance<
‎      typename boost::mpl::begin<FieldSequence>::type,
‎      typename boost::mpl::find<
‎        FieldSequence,
‎        field<key_at_pos::value>
‎      >::type
‎    >::value;
‎  return t.template get<M>();
‎}

The code is simpler than it might appear at first glance: if the key at the position Pos of CompositeKey is of the form key<n>, M is simply the distance from the beginning of FieldSequence to the type key<n>. The expression key_at_pos::value takes advantage of the fact that each key<n> is derived from boost::mpl::int_<n>.

The different overloads of find are entirely similar to the one we have just seen. Boost.Preprocessor allows us to generate the overloads without code repetition:

‎#define FIND_PARAM_MACRO(z,n,data) \
‎const field<BOOST_PP_CAT(N,n)>& BOOST_PP_CAT(f,n)
‎
‎#define FIND_FIELD_MACRO(z,n,data) field<BOOST_PP_CAT(N,n)>
‎
‎#define FIND_VALUE_MACRO(z,n,data) BOOST_PP_CAT(f,n).x
‎
‎#define FIND_GET_MACRO(z,n,data) \
‎get<composite_key_type,fields_type,boost::mpl::int_<n> >(t)
‎
‎#define DEFINE_FIND(num_fields) \
‎template<BOOST_PP_ENUM_PARAMS(num_fields,int N)> \
‎const element* find( \
‎  const container& c, \
‎  BOOST_PP_ENUM(num_fields,FIND_PARAM_MACRO,~)) \
‎{ \
‎  typedef cover< \
‎    cover_number< \
‎      boost::mpl::vector_c<int, \
‎        BOOST_PP_ENUM_PARAMS(num_fields,N) \
‎      > \
‎    >::value \
‎  >                                                tag; \
‎  typedef typename container::index<tag>::type     index_type; \
‎  typedef typename index_type::key_from_value      composite_key_type; \
‎  typedef boost::mpl::vector< \
‎    BOOST_PP_ENUM(num_fields,FIND_FIELD_MACRO,~) \
‎  >                                                fields_type; \
‎ \
‎  const index_type& i=c.get<tag>(); \
‎  boost::tuple< \
‎    BOOST_PP_ENUM_PARAMS(num_fields,int BOOST_PP_INTERCEPT) \
‎  > t=boost::make_tuple( \
‎    BOOST_PP_ENUM(num_fields,FIND_VALUE_MACRO,~) \
‎  ); \
‎  typename index_type::iterator it=i.find( \
‎    boost::make_tuple( \
‎      BOOST_PP_ENUM(num_fields,FIND_GET_MACRO,~) \
‎    ) \
‎  ); \
‎  if(it!=i.end())return &*it; \
‎  else return 0; \
‎}
‎
‎DEFINE_FIND(1)
‎DEFINE_FIND(2)
‎DEFINE_FIND(3)
‎DEFINE_FIND(4)

A complete implementation program is provided.

Although we have metaprogrammed the code for selecting and using the indices of an n-attribute container, still the definition of container was done manually. The following is a very tough challenge for the reader: program a metafunction that accepts an integer n and produces a suitable multi_index_container instantiation providing full query coverage for a struct with attributes x1,...,xn. Note that solving this challenge implies metaprogramming the permutation cover algorithm we devised some entries ago.

Filling up at dawn

2008-11-04T01:01:00.002+01:00

I recently read in a free daily paper that it is advisable to fill up the car tank early in the morning rather than at noon because the gasoline will be usually cooler at that time of the day and thus denser, so we are pumping more grams of gas per liter paid. Does this really make a difference?

The volumetric thermal expansion coefficient of gasoline is 950·10^-6/°C, which means that when the gasoline is colder we obtain an excess mass of

(m_cold − m_hot)/m_hot = 950·10^-6ΔT.

For a rather large ΔT of 20 °C the gain would be around 2%. Take into account, however, that station tanks are buried a couple of meters underground, which damps the temperature oscillations experienced at the surface. This damping effect can be estimated as:

ΔT(z) = ΔT_surface·e^−z/D,
D = √(2K/ω),

where K is the soil thermal diffusivity, with typical values for natural soils ranging from 0.2 to 0.8 10^-6m²/s. Setting z = 2 m and ω = 2π rad/day, we have

1.9·10^-12ΔT_surface ≤ ΔT_{station tank} ≤ 1.4·10^-6ΔT_surface,

that is, the temperature of the station tank is basically constant across the day. So, filling up the car tank in the cool hours of the day does not seem to make any difference.

There is a measurable gain, though, when we compare car tank filling during different seasons of the year. Here, the temperature oscillations are higher and, most importantly, the damping effect of the soil is much lesser. Taking

ΔT_surface = 30 °C,
z = 2 m,
ω = 2π rad/year,

we have

7.3 °C ≤ ΔT_{station tank} ≤ 14.8 °C,

which translates to

0.7% ≤ (m_cold − m_hot)/m_hot _t ≤ 1.4%.

The potential impact as a cost-saving measure is, in any case, very unspectacular.

C preprocessor tail recursion

2008-10-31T07:35:00.001+01:00

The C preprocessor gives the false impression of being functional in nature, but its capacities are more limited than expected; for instance, function reentrancy does not work: consider the following attempt at recursively defining the factorial function using Boost.Preprocessor:

‎#define FACTORIAL(n) \
‎BOOST_PP_IF( \
‎  n, \
‎  BOOST_PP_MUL(n,FACTORIAL(BOOST_PP_DEC(n))), \
‎  1)

The expression FACTORIAL(5) does not expand to 120, however, but to

‎BOOST_PP_TUPLE_ELEM_3_0 BOOST_PP_IIF_BOOST_PP_BOOL_FACTORIAL(4)(
‎BOOST_PP_WHILE_2, (0, 5, FACTORIAL(4)) BOOST_PP_TUPLE_EAT_3)(BOO
‎ST_PP_MUL_P, BOOST_PP_MUL_O, (5, 5, BOOST_PP_DEC_FACTORIAL(4)))

The rules of the C preprocessor mandate that during the expansion of the function macro FACTORIAL internal occurrences of this token are not further expanded, thus ruling out the very core of recursion. A simpler example illustrates this:

‎#define A B
‎#define B A
‎
‎A

The last line of the snippet above does not enter into infinite recursion, but simply causes the expansion rules A→B→A to be applied, stopping here because A is met inside the expansion of A itself.

Despite these serious limitations, Boost.Preprocessor is able to provide constructs for bounded iteration implemented with very sophisticated techniques. Let us see how to use BOOST_PP_WHILE to simulate tail recursion, a special form of recursion in which recursivity is applied exactly at the end of the recursive function body. The usual formulation of the factorial function is not tail recursive:

‎int factorial(int n)
‎{
‎  if(n>0)return n*factorial(n-1);
‎  else   return 1;
‎}

But it can be rewritten to be so as follows:

‎int factorial(int n,int res=1)
‎{
‎  if(n>0)return factorial(n-1,n*res);
‎  else   return res;
}

The idea is that the recursive function passes all the relevant contextual information as part of the recursive call, thus delegating the work that otherwise should be done after the recursion returns; in this way, recursing (except when terminating recursion) is exactly the last thing that the function does. Let us know apply the following trasformation f → f' to a tail recursive function f:

If f stops and returns res, f' returns the sequence (stop,res).
If f recursively calls f(args), f' returns the sequence (iterate,args).

That is, f' does not recurse but instead returns the computation state if recursion is to occur, or the final return value otherwise. The following pseudocode illustrates the transformation process for factorial:

‎factorial_body(int n,int res=1)
‎{
‎  if(n>0)return (iterate,n-1,n*res);
‎  else   return (stop,res);
‎}

f' is no longer recursive, but it can be used to compute f as follows:

‎recurse(body,args)
‎{
‎  state=(iterate,args);
‎  while(state[0]!=stop){
‎    state=body(state.args);
‎  }
‎}

The interesting fact about this process is that recursion (and so reentrancy) has been eliminated in favor of iteration: f' is called repeatedly until it signals termination. We can very easily implement this idea using Boost.Preprocessor:

‎#define RECURSE(body,args) \
‎RECURSE_RES(BOOST_PP_WHILE(RECURSE_PRED,RECURSE_OP,(body)ITERATE args))
‎
‎#define RECURSE_PRED(d,state) \
‎BOOST_PP_SEQ_ELEM(1,state)
‎
‎#define RECURSE_OP(d,state) \
‎(RECURSE_BODY(state))RECURSE_BODY(state)(BOOST_PP_SEQ_REST_N(2,state))
‎
‎#define RECURSE_BODY(state) \
‎BOOST_PP_SEQ_ELEM(0,state)
‎
‎#define RECURSE_RES(state) \
‎BOOST_PP_SEQ_ELEM(2,state)
‎
‎#define STOP (0)
‎
‎#define ITERATE (1)

which allows us to implement our factorial macro in a reasonably straightforward manner:

‎‎#define FACTORIAL(n) RECURSE(FACTORIAL_BODY,(n)(1))
‎
‎#define FACTORIAL_BODY(args) \
‎FACTORIAL_BODY_IMPL( \
‎  BOOST_PP_SEQ_ELEM(0,args),BOOST_PP_SEQ_ELEM(1,args))
‎
‎#define FACTORIAL_BODY_IMPL(n,res) \
‎BOOST_PP_IF( \
‎  n, \
‎  ITERATE(BOOST_PP_DEC(n))(BOOST_PP_MUL(n,res)), \
‎  STOP(res))

Arguments are packed using Boost.Preprocessor sequences. A complete implementation program is provided. Unfortunately, the process seems to be very onerous for common implementations of the C preprocessor, as internal limits are generally reached with arguments greater than 5. This problem could deserve further investigation in a later entry.

Law-abiding miracles

2008-10-27T03:58:00.001+01:00

Newton believed that, if left to the action of gravitational forces alone, the Solar System must be unstable, and only the continuous action of God, counteracting the small perturbations here and there, keeps the system together and wonderfully regular:

The six primary planets are revolved about the sun in circles concentric with the sun, and with motions directed towards the same parts, and almost in the same plane. Ten moons are revolved about the earth, Jupiter and Saturn, in circles concentric with them, with the same direction of motion, and nearly in the planes of the orbits of those planets; but it is not to be conceived that mere mechanical causes could give birth to so many regular motions[...] This most beautiful system of the sun, planets, and comets could only proceed from the counsel and dominion of an intelligent and powerful Being.

This curious picture of God as a careful supervisor of heavenly motions would only be compatible with the same physical laws that Newton crafted if God Himself would possess a physical body with mass and energy (and large quantities of them, for that matter): classical (and relativistic) mechanics is deterministic and does not allow for any degree of freedom that a deity could take advantage of to direct the evolution of a physical system one way or another. In this sense, the action of God as a stabilizer of the Solar System is a sustained miracle, a physical phenomenon incompatible with the laws of physics. This very definition of miracle is logically contradictory to begin with, at least given the premise of a deterministic universe.

On the other hand, if we consider non-deterministic physical theories, there is room for a concept that we could aptly name "law-abiding miracle": for instance, quantum mechanics tunneling allows for the occurence of phenomena that we generally deem impossible, such as wall crossing and other ghostly effects, albeit the probability of their occurrence at a macroscopic level is exceedingly small. So small in fact that it is not reasonable to expect that a macroscopic tunneling event takes place during the entire lifetime of the Universe, let alone in the presence of humans. But unlikely is not impossible: a quantum god can in principle miraculously interfere with the business of the physical world without violating the rules that govern it.

There is, nonetheless, some degree of circularity with this notion of law-abiding miracles within a probabilistic physical theory. If the frequency of such miracles is too high (as it would be the case with the scenario envisioned by Newton of a solar system being constantly tweaked to keep it stable), such divine interventions would show up in the statistical calculations made by physicists when experimentally testing the theory, rendering results highly inconsistent with the theoretical predictions; and this can in the end falsify the theory. The insistence of a deity to meddle in the normal course of events would block humans from acquiring knowledge of the laws of Nature.

Circles, ellipses, assistants and employees

2008-10-23T22:30:00.002+02:00

The well-known Circle-ellipse Problem states that, in the world of OO, a circle is not an ellipse because both entities do not behave the same: for instance, one can arbitrarily stretch an ellipse, but stretching a circle would break the essential property of this shape, namely that all its points are at the same distance from the center. On the other hand, other subtype relationships perfectly match the real-life scenarios they model: for instance, an assistant is certainly a kind of employee, and so an associated Assistant type can be modeled as a subtype of Employee without any problem whatsoever. What is the difference between these two modelizations that make the latter succeed where the former fails?

We introduce notation to talk about some concepts around the notion of "object" in OO. An object x has a definite lifetime T(x) = [t₀,t₁] during which the internal state of the object might (or might not) change according to the operations performed on it. We denote by S(x,t) the state of x at time t. Assume that we have a type Ellipse modelling the mathematical concept of ellipse, which we can identify by the set E containing all the ellipses in R²: we intuitively expect from Ellipse that the states of objects of this type can be univocally mapped to elements of E:

for all x of type Ellipse, t in T(x), there is an ellipse e in E such that e ~ S(x,t),

or, put more succintly, we can say that

the states of Ellipse objects are ellipses (modulo isomorphism).

Similarly, the states of Circle objects are circles. A Circle cannot be a subtype of Ellipse because, according to Liskov Principle, a Circle does not behave as an Ellipse: there are mutable operations on Ellipse (like the streching example mentioned at the beginning) such that the final state of the modifed object can never be a circle, so we cannot provide an equivalent operation for Circle.

Now, suppose we have types Employee and Assistant: here the subtyping relationship works because all the state changes of an employee (salary raise, relocation, etc.) are compatible with the state changes of an assistant. And this points right to the crux of the matter: the real-life "employee" concept refer to a mutable entity whose state can change as part of the interactions of the employee with the rest of the company; Employee models this concept by making the object state changes reflect the changes of the modelled employee:

the states of Employee objects are employee states.

Confront this with

the states of Ellipse objects are ellipses.

Ellipse states are ellipses, not ellipse states (ellipses do not have any mutable state). This is why the the subtypying relationship Ellipse/Circle does not work while Employee/Assistant does. In general, for a concept T < T' relationship to be representable by types T < T', the states of T and T' objects must model states of T and T' entities, respectively.

Multicolumn querying

2008-10-20T00:08:00.002+02:00

Suppose we have a database table equipped with a multicolumn (also called concatenated or complex) index like this:

‎CREATE TABLE user (last_name varchar, first_name varchar)
‎CREATE INDEX user_index ON user (last_name, first_name)

Most platforms take advantage of user_index in queries involving both user attributes or just last_name:

‎SELECT * from user WHERE first_name="Robert" AND last_name="Grant";
‎SELECT * from user WHERE last_name="Smith";

while, on the other hand, queries like the following:

‎SELECT * from user WHERE first_name="Mary";

trigger a full table scan or use the index in a suboptimal fashion. The simple reason behind this is that database indexes are usually implemented with B-trees or similar structures that sort the data by the values of the indexed column: for a multicolumn index on (c₁,...,c_n), sorting is performed according to the lexicographical order associated to the n columns, with c_i taking precedence over c_j if i < j. So, in our example the index sorts the rows by last_name and, within clusters with equal last_name, data is further sorted by first_name. This explains why querying on last_name benefits from the index (rows are sorted by the first attribute to begin with) while querying on first_name does not. Database designers learn to carefully consider the column specification order when creating multicolumn indexes.

So, in general, an index on (c₁,...,c_n) serves to accelerate queries of the form

c₁ = v₁ AND c₂ = v₂ AND ··· AND c_i = v_i

for all i ≤ n. The question arises: how many multicolumn indexes, and how, must be created so that all possible attribute combinations are covered by some index? For instance, if we have three attributes:

‎CREATE TABLE user (last_name varchar, first_name varchar, age integer)

three indexes suffice:

‎CREATE INDEX user_index1 ON user (last_name, first_name, age)
‎CREATE INDEX user_index2 ON user (first_name, last_name, age)
‎CREATE INDEX user_index3 ON user (age, last_name, first_name)

(You might find it entertaining to check out manually that all the 7 different attribute combinations are covered.) It turns out that for n attributes, which generate 2ⁿ − 1 combinations, the minimum number of indexes needed is

C(n,floor(n/2)) = n!/(ceil(n/2)!·floor(n/2)!),

where C(a,b) denotes the binomial coefficient a choose b. The values of C(n,floor(n/2)) for n = 1,2,... are: 1, 2, 3, 6, 10, 20, 35, 70, 126,... Using Stirling's approximation, the sequence can be proved to be O(2ⁿ/√n). The mathematical formulation of this problem has been studied in a previous entry. Also, an algorithm for generating the index set is provided.

On necessity in McGinn

2008-10-16T05:58:00.005+02:00

Colin McGinn challenges in his book Logical Properties the standard treatment of modality and proposes some rather provocative and unorthodox ideas on the nature of necessity and associated modal notions. There is much to argue about McGinn's positions, but I would like to focus on three particular issues:

The main objection of McGinn to the possible world semantics is that this approach reduces modal notions to quantifications over universes of possible worlds, and whether a world is possible or not is already dependent on the modal concept of possibility, thus engendering a circular reasoning. As I see it, the flaw in this argument stems from the wrong assumption that characterizing a world W as possible is a modal business; although McGinn does not put it explicitly, he seems to think that the statement "W is a possible world" is equivalent to something like "possibly W exists". But actually there is nothing modal about the concept of "possible world": a possible world is just a consistent (i.e. logically non-contradictory) state of affairs, of which the actual world we live in is just an example. Rephrase the whole possible worlds semantics formulation replacing "possible world" with "consistent world" and all the apparent circularities vanish away.

Later on, McGinn proposes a copula modifier theory to explain modality: basically, this theory contends that modal notions qualify the way that objects hold properties. So, "Socrates must be a man" means that Socrates holds the property of being a man in the mode of necessity. McGinn immediately realizes that adscribing modality to the copula of statements cannot cover uses of modality qualifying entire statements, as in "it is impossible that 2 + 2 = 5" (which, incidentally, is how modality operators work in formal modal logic). He solves the problem by postulating that in these cases the missing copula to qualify is that adscribing the truth property to the proposition at hand: so "it is impossible that 2 + 2 = 5" can be rephrased as "the proposition '2 + 2 = 5' cannot be true", which, according to the copula modifier theory, means that "2 + 2 = 5" holds the property of being true in the mode of impossibility. In the process, McGinn has forced us to accept propositions as first-class entities to talk about, an ontological extension that some (notably Quine) reject commiting to. Even admitting this unexpected guest, it is not clear to me how this theory sheds any light on the nature of modality: one can as well go the reverse way and formulate predicate-qualifying modal sentences as sentence-level modal utterances: for instance, "Socrates must be a man" is just a way of saying "necessarily Socrates is a man", which process does not involve the introduction of propositions and the "truth predicate". Why McGinn prefers one rewriting rule to the converse is unknown to me --furthermore, McGinn's formulation poses some perplexities when iterated modality is considered, as he readily acknowledges.

Finally, having stated that modality affects the way objects hold properties, McGinn renounces to further analyze what the exact nature of this relationship between objects and properties is. He gives in passing a Tarski-style semantic formulation of modality along the lines

x satisfies 'is necessarily F' iff x necessarily satisfies 'F',

only to admit that this does not add anything to our understanding of modality. Modal truth is left as a sort of primitive, epiphenomenal concept for which even a special entity status is set up within the confines of McGinn's exuberant ontology. It is difficult to see how this situation represents an improvement over the standard modal theory.

Points inside a polygon and the Residue Theorem

2008-10-12T07:51:00.001+02:00

A usual way to determine whether a 2D point p lies inside a simple polygon P consists in tracing a ray from p to infinity and computing how many times this ray intersects P: p is inside P iff the number of intersections is odd. We investigate a different technique for solving the problem based on complex analysis, only to finally see that this new method is really equivalent to the classical ray algorithm.

We consider p and P to be represented in the complex plane and assume, without loss of generality, that p = 0. The Residue Theorem tells us that

∫_P (1/z) dz

is null if the point 0 (the only pole of 1/z) is outside P, and

2πi Res_z=0(1/z) = 2πi

if 0 lies inside P. So we need only evaluate the integral ∫_P (1/z) dz; as P is a loop sequence of segments going counterclockwise like this

P = (p₁,p₂), (p₂,p₃), ... , (p_N−1,p_N), (p_N,p₁),

and ln z is an antiderivative of 1/z, we have that

∫_P (1/z) dz = ∑ ∫_{(p_i,p_i+1)} (1/z) dz =
= ∑ (ln p_i+1 − ln p_i).

But now, this sum seems to be unconditionally null:

∑ (ln p_i+1 − ln p_i) =
= ln p₂ − ln p₁ + ln p₃ − ln p₂ + ln p₄ − ln p₃ + ... + ln p₁ − ln p_N = 0.

Each term is cancelled out by some other with different sign, so the total sum is, or seems to be, zero! Obviously we have made some mistake during the process, and in fact the wrong assumption is this: the equality

∫_{(p_i,p_i+1)} (1/z) dz = ln p_i+1 − ln p_i

only holds if the derivative of ln z is 1/z along all the segment (p_i,p_i+1) (more precisely, at a region containing the segment): but ln z, or, strictly speaking, the principal branch of ln z, has a cut at the negative real axis, so the equality above fails to hold when (p_i,p_i+1) intersects this cut. When this is the case, we can evaluate the integral by using a different branch of the complex logarithm so that (p_i,p_i+1) does not intersect the new cut. It is not hard to see that doing so yields

∫_{(p_i,p_i+1)} (1/z) dz = ln p_i+1 − ln p_i + σ2πi,

where ln denotes again the principal branch of the complex logarithm and σ is +1 or −1 depending on whether (p_i,p_i+1) crosses the negative real axis southwards or northwards. Having the situation with the negative real axis covered, we can easily deduce that

∫_P (1/z) dz =(s − n)2πi,

where s is the number of segments of P crossing the negative real axis southwards, and n corresponds to those crossing northwards. And, in the end, this is but a simple reformulation of the old ray-based algorithm.

Generating permutation covers: part II

2008-10-08T04:17:00.001+02:00

In a prior entry, we reduced the problem of generating a minimal permutation cover of a set S of N elements to that of constructing a minimal tuple cover of S₀ U ... U S_M, where S_i is the set of all subsets of S with cardinality i, M = (N − 1)/2 and N is odd. We do this construction recursively on m = 0,...,M.

Suppose then that we have a minimal tuple cover T of S₀ U ... U S_m, i.e. a minimal set of m-tuples of different elements jointly covering all subsets of S with cardinality ≤ m, and we want to extend T to a minimal set T' of (m + 1)-tuples jointly covering S₀ U ... U S_m U S_m_{+ 1}. We state without proof the following

Lemma. If m ≤ N/2, a miminal tuple cover of S₀ U ... U S_m has cardinality C(N,m).

So, |T| = C(N,m) and we want to extend T to a set T' by assigning to each tuple τ of T one or more tuples of the form τa, a in S − range(τ), in such a way that all subsets of S_m_{+ 1} are covered and the final number of elements of the cover |T'| is C(N,m + 1).

When extending a tuple τ of T we will disregard the order of its elements, so basically we are identifying τ with range(τ); this reduction still yields C(N,m) different elements, since every element of S_m must be covered by some tuple of T. So the extension mapping T → T' relates elements of S_m to elements of S_m_{+ 1}, and in fact can be regarded as a subset of the graph induced by inclusion between elements of these two sets.

It is easy to see that each element of S_m is the source of (N − m) arrows, whereas m + 1 arrows arrive at each element of S_m_{+ 1}. Our job is then to select C(N,m + 1) arrows from the diagram above so that all the elements from the source and the destination set are met by at least one arrow. This selection process has to maintain some balance so that no source or destination element is left unattended: for instance, if we select for each element of S_m_{+ 1} an arbitrary arrow arriving at it there is the possibility that some elements of S_m are not visited by any of the selected arrows. The following is a balanced selection criterion: given a fixed injection f : S_m → S_m_{+ 1} and function g : S_m_{+ 1} → S_m, both compatible with the inclusion relationship between S_m and S_m_{+ 1} (that is, X is a subset of f(X) and Y a superset of g(Y)), we select an arrow X → Y iff

Y = f(X) or (X = g(Y) and Y is not in f(S_m)).

It is obvious that this criterion does not leave any X or Y unvisited. The number of selected arrows coincides with the number of elements in S_m_{+ 1}, which is C(N,m + 1) as required. In the following we suppose, without loss of generality, that S is the set {0,...,N − 1}. Constructively defining an inclusion compatible injection from S_m to S_m_{+ 1} is not a trivial task, but fortunately for us a paper from Pudlák, Turzík and Poljak provides the definition of such an injection f along with an algorithm χ : S_m_{+ 1} → {true,false} that checks whether a given Y belongs to f(S_m). We adopt the following definition for g:

g(Y) := Y − max(Y),

which leads to this algorithm for generating Τ' from Τ:

extend(N,Τ)
Τ' ← Ø
for every τ in Τ
····X ← range(τ)
····a ← f(X) − X {the one element added to X by f}
····Τ' ← Τ' U {τ·a}
····for i = max(X) + 1 to N − 1
········Y ← X U {i}
········if not χ(Y) then
············Τ' ← Τ' U {τ·i}
········end if
····end for
end for

Note that the double loop over (X,i) is designed in such a way that it only visits the X → Y arrows where X = g(Y), which is maximally efficient and saves us the need to explicity compute g. The complete Τ₀ = Ø → Τ₁ → ··· → Τ_M process can be inlined to avoid generating the intermediate Τ_i covers.

tuple-cover(m,N,τ,Τ) {initial call: tuple-cover(0,N,Ø,Τ) with Τ empty}
if m = (N − 1)/2 then
····Τ ← Τ U {τ}
else
····X ← range(τ)
····a ← f(X) − X
····tuple-cover(m + 1,N,τ·a,Τ)
····for i = max(X) + 1 to N − 1
········Y ← X U {i}
········if not χ(Y) then
············cover(m + 1,N,τ·i,Τ)
········end if
····end for
end if

In order to leverage the tuple cover algorithm to construct a minimal permutation cover, the only missing piece is finding a bijection d : S_M → S_M such that X ∩ d(X) = Ø for all X in S_M, as we already saw. The aforementioned paper from Pudlák et al. provides also such a function (which, in fact, is used for the construction of f). Having all the necessary components, the following is the full algorithm for constructing a minimal permutation cover on S = {0,...,N − 1}.

cover(N)
if N is even then
····N' ← N − 1
else
····N' ← N
end if
Σ ← Ø
Τ ← Ø
tuple-cover(0,N',Ø,Τ)
for every τ in Τ
····find τ' in Τ with range(τ') = d(range(τ))
····a ← S − range(τ) − range(τ')
····σ ← τ·a·reverse(τ')
····if N is even then
········Σ ← Σ U {N'·σ}
········Σ ← Σ U {σ·N'}
···· else
········Σ ← Σ U {σ}
····end if
end for

A C++ implementation of the algorithm is available (Boost used). The following are the different covers generated for |S| = 1,...,8.

\|S\|	permutation cover
1	(a)
2	(ab), (ba)
3	(acb), (bac), (cba)
4	(acbd),(dacb),(bacd),(dbac),(cbad),(dcba)
5	(acedb),(baedc),(bdaec),(bedca),(cbaed),(cdbae),(cebad),(daceb),(decab),(eadbc)
6	(acedbf),(facedb),(baedcf),(fbaedc),(bdaecf),(fbdaec),(bedcaf),(fbedca),(cbaedf),(fcbaed), (cdbaef),(fcdbae),(cebadf),(fcebad),(dacebf),(fdaceb),(decabf),(fdecab),(eadbcf),(feadbc)
7	(acegfdb),(baegfdc),(bdagfec),(bdfagec),(bdgfeca),(bedagfc),(befdcag), (begdcaf),(bfaegdc),(bfgecad),(bgafced),(cbagfed),(cdbagfe),(cdfbage), (cdgbafe),(cebagfd),(cefbagd),(cegbafd),(cfbaged),(cfgeadb),(cgbfdae), (dacgfeb),(decbagf),(defcagb),(degcafb),(dfacgeb),(dfgceab),(dgafbec), (eadcgfb),(efadbgc),(efgdabc),(egadbfc),(facegdb),(fgaebdc),(gacfdeb)
8	(acegfdbh),(hacegfdb),(baegfdch),(hbaegfdc),(bdagfech),(hbdagfec),(bdfagech), (hbdfagec),(bdgfecah),(hbdgfeca),(bedagfch),(hbedagfc),(befdcagh),(hbefdcag), (begdcafh),(hbegdcaf),(bfaegdch),(hbfaegdc),(bfgecadh),(hbfgecad),(bgafcedh), (hbgafced),(cbagfedh),(hcbagfed),(cdbagfeh),(hcdbagfe),(cdfbageh),(hcdfbage), (cdgbafeh),(hcdgbafe),(cebagfdh),(hcebagfd),(cefbagdh),(hcefbagd),(cegbafdh), (hcegbafd),(cfbagedh),(hcfbaged),(cfgeadbh),(hcfgeadb),(cgbfdaeh),(hcgbfdae), (dacgfebh),(hdacgfeb),(decbagfh),(hdecbagf),(defcagbh),(hdefcagb),(degcafbh), (hdegcafb),(dfacgebh),(hdfacgeb),(dfgceabh),(hdfgceab),(dgafbech),(hdgafbec), (eadcgfbh),(headcgfb),(efadbgch),(hefadbgc),(efgdabch),(hefgdabc),(egadbfch), (hegadbfc),(facegdbh),(hfacegdb),(fgaebdch),(hfgaebdc),(gacfdebh),(hgacfdeb)

In a later entry we will see a practical application of permutation covers in the context of database querying.

User attention sharing

2008-10-04T21:22:00.004+02:00

CPU is a precious resource that must be shared between the processes concurrently running on a machine at a given time. Some operating systems of the past implemented this sharing using a technique called collaborative multitasking, which is based on the processes' explicitly relinquishing CPU control to the OS so that this control be passed on to other tasks. This can be very inefficient, as sharing cannot be controlled by the operating system, but depends on the good judgement of third party application programmers; furthermore, a faulty process that fails to yield CPU or hangs will collapse the whole system. So, modern OSs implement preemptive multitasking, in which the system scheduler can give or take CPU control without any special help from the processes. Usual implementations of preemptive multitasking group processes into two categories: IO-bound processes, which consume most of the time waiting for user input or dealing with secondary storage, and CPU-bound processes, which take up as much computing power as they are given. Some scheduling algorithms reward the least CPU-hungry processes by giving them higher priority to regain CPU.

In a desktop machine, user attention is as much of a precious shared resource as CPU is: concurrently living applications try to attract the user focus by displaying some visual effects or directly putting themselves on the forefront. In the early, more naïve days, applications could freely jump into the foreground of their own accord, which could quickly degrade into an all-out fight (specially if malware was involved) for the user's attention, or "attention thrashing", and render the computer unusable. This situation is reminiscent of the first candid collaborative multitasking systems; incidentally, attention thrashing led some OS vendors to change the behavior of some of their APIs to try to restrict the applications' freedom to steal focus from each other, though workarounds abound. In the web browser arena, attention thrashing can be even more daunting, and most users end up entirely blocking pop-ups and other forms of focus stealing.

As I see it, user attention should then be explicitly taken care of by the operating system in the same way as CPU control is; in fact, we can learn from the history of multitasking and move from the current "collaborative" situation to system-controlled scenarios where applications request attention and are given it based on attention availability and some measure of application fairness. So, applications that request attention very rarely and do so for justified motives (where this last attribute can be estimated indirectly by, say, recording the amount of time spent by the user on the application when this pops up) are ranked higher and can be given more leeway to call on the user. On the other hand, noisy applications are rapidly ranked low and will not be able to steal focus unless the user explicitly clicks on them.

Mass exhalation

2008-09-30T01:05:00.002+02:00

How much bodily mass do we give out through breathing? The chemical reaction for aerobic respiration is

C₆H₁₂O₆ + 6O₂ → 6CO₂ + 6H₂O

The amount of carbon dioxide we exhale is then

6 mol CO₂ / mol glucose = 1.47 g CO₂ / g glucose.

Processing of fatty acids yields a different amount of exhaled CO₂. For instance, the oxidation of a saturated fatty acid is expressed as

C_nH_2nO₂ + (3n/2 − 1) O₂ → n CO₂ + n H₂O

yielding

44n/(14n + 32) g CO₂ / g C_nH_2nO₂.

This last expression also serves approximately for unsaturated fatty acids, since these are missing only a few lightweight hydrogen atoms with respect to the saturated acid with the same amount of carbon. Assuming most consumed fat is in the form of palmitate, we have a carbon dioxide yield of

2.75 g CO₂ / g fat;

Carbohydrates produce 4 kcal per gram, whereas fatty acid processing yields 9 kcal per gram of fat. While at rest, the human body takes 30% of the energy from carbohydrates and 70% from fat; so, this results in a combined production of carbon dioxide of

0.324 g CO₂ / kcal.

For a daily energy consumption rate around 2,500 kcal, 810 grams of carbon dioxide are then expelled through breathing. There is a experimental study by the USDA on the subject that unfortunately I haven't been able to find, but some sources cite it as reporting a value beween 445 and 450 liters of CO₂ (881-891 g) per day, which is in good agreement with our calculations.

The value calculated is not the answer to our original question about exhaled bodily mass: much of the oxygen in the CO₂ produced comes from the air rather than catabolized substances, specially in the case of fat. We will assume that all the oxygen coming from the burned substance goes to CO₂: under these conditions, the amount of bodily mass expelled through breathing is

0.933 g / g glucose,
0.875 g / g fat,

or, using the assumed ratios of carbohydrate and fat consumption,

0.138 g / kcal,

which yields 345 g for a typical daily energy consumption of 2,500 calories.

Generating permutation covers: part I

2008-09-26T05:50:00.003+02:00

We investigate ways of efficiently generating minimal permutation covers for a set S of N elements. As we have already proved, such covers have size C(N,floor(N/2)). The following intermediate results reduce the problem to simpler terms.

Theorem. Assume that N is even and greater than zero. Let a be an an arbitrary element of S and Σ a minimal permutation cover of S − {a}. The set Σ' defined as

Σ' = {aσ, σa : σ in Σ}

is a minimal permutation cover of S.

Proof. Σ' is a permutation cover: let X be a subset of S; if a is not in X, X is covered by some permutation σ of Σ and a fortiori by σa in Σ'; if a belongs to X, some σ in Σ covers X − {a} and then aσ in Σ' covers X. Σ' is minimal as |Σ'| = 2|Σ| = 2 C(N − 1, floor((N − 1)/2)) = 2 C(N − 1, N/2 − 1) = C(N,N/2) = C(N,floor(N/2)), where the equalities depend on basic properties of binomial coefficients and the fact that N is even.

This result reduces the problem to N odd. The Hasse diagram of the powerset of S when N is odd looks like this (the particular case depicted is N = 5):

We name the layers of the diagram S₀,...,S_N. S_i consists of the subsets of S with cardinality i. We have |S_i| = C(N,i). A permutation covers exactly one element of each S_i. Similarly, a tuple τ of n non-repeating elements covers exactly one element of each of S₀,...,S_n. So, a minimal tuple cover of S₀,...,S_n can be proved to have as many elements as the largest layer covered. The next result shows how to extend a minimal tuple cover of S₀,...,S_{(N − 1)/2}, which has cardinality C(N,(N − 1)/2) = C(N,floor(N/2)) and spreads over the bottom half of the powerset of S, to a minimal permutation cover of S.

Theorem. Assume that N is odd and set M = (N − 1)/2. Let Τ be a minimal set of M-tuples of non-repeating elements covering S₀ U ... U S_M, and d : S_M → S_M a bijection with X ∩ d(X) = Ø for all X in S_M. The set Σ defined as

Σ = { τ·a·reverse(τ') : τ,τ' belong to Τ, range(τ') = d(range(τ)), a = S − range(τ) − range(τ')}

is a minimal permutation cover of S.

Proof. We see first that to each τ in Τ there corresponds exactly one permutation σ_τ of the form τ·a·reverse(τ') described above: since Τ is minimal each element of S_M is covered by one and only one element of Τ; so, for each τ there is exactly one τ' covering d(range(τ)). We prove now that Σ covers S: let X be a subset of S; if |X| ≤ M, X is covered by some τ in Τ and hence also by the associated σ_τ; if |X| > M, the set Y = S − X has |Y| ≤ M and is thus covered by some τ' in Τ; as we have seen, τ' univocally determines a permutation σ = τ·a·reverse(τ'), and σ obviously covers X, as all the elements of Y are packed at the tail of the permutation. Finally, Σ is minimal, since |Σ|= |Τ| and we saw before that |Τ| = C(N,floor(N/2)).

The mapping τ → σ_τ just defined can be described in a more algorithmic fashion as follows:

Find the τ' in Τ covering d(range(τ)).
Let a be the only element of S not present in τ nor τ'.
σ_τ := τ·a·reverse(τ').

Note that the last theorem extends easily to N even: in this case, the permutations have the form τ·reverse(τ'), i.e. there is no middle element a. However, it is more efficient to reduce the case N even to N − 1 using the first theorem of this entry.

Now we are left with the tasks of designing an algorithm to generate a tuple set Τ covering S₀ U ... U S_M and finding a suitable bijection d : S_M → S_M.

Profiling stable_vector

2008-09-22T02:00:00.006+02:00

The two random access containers of the C++ standard library, std::vector and std::deque, have the same big-O complexity as stable_vector for the following operations:

push_back: amortized constant
middle insert: O(n)
operator[]: constant

Yet the differences in actual performance vary wildly. I have written a profile program that measures the time taken by the three containers at these operations for two different value types: int and a "heavy" user defined type bigobj holding an int and an std::string. The program was built with release mode settings with MSVC++ 8.0 and GNU GCC 3.4.4 for Cygwin/MinGW. Resulting times are summarized in the following table: values are normalized to the minimum among the three containers.

MSVC++ 8.0
	`int`			`bigobj`
	`std:: vector`	`std:: deque`	`stable_ vector`	`std:: vector`	`std:: deque`	`stable_ vector`
`push_back`	1.00	4.06	21.37	1.00	1.09	1.26
`insert`	1.00	4.24	5.29	1.50	1.73	1.00
`operator[]`	1.00	3.35	1.01	1.00	3.27	1.09
GNU GCC 3.4.4
	`int`			`bigobj`
	`std:: vector`	`std:: deque`	`stable_ vector`	`std:: vector`	`std:: deque`	`stable_ vector`
`push_back`	1.09	1.00	91.57	1.08	1.00	1.14
`insert`	1.00	1.30	5.17	21.38	21.36	1.00
`operator[]`	1.00	41.04	1.11	1.00	81.68	2.03

For small objects like int, std::vector is, as expected, the fastest or very close to the fastest at every operation. push_back takes much longer for stable_vector than std::deque presumably because the latter allocates element chunks thus amortizing allocation time whereas stable_vector performs an allocation for every new node. The poor performance of stable_vector at insert can be attributed to the operation of adjustment of "up" pointers with each insertion operation, which is considerably slower than the simple int copying incurred by the other two containers. Finally, note that random access is very slow for std::deque because it involves some non-trivial arithmetical calculations, whereas std::vector and stable_vector only perform simple pointer offsets and indirections.

When storing heavy objects with expensive copy, the situation changes somewhat. The fact that stable_vector does not internally copy elements upon resizing or shrinking improves its relative performance with respect to the other two containers at push_back, and specially at insert, where stable_vector is the fastest.

This basic analysis suggests that the main performance handicap of stable_vector stems from the fact that it is a node-based container, as opposed to the chunk-based nature of std::vector and std::deque. This drawback might be alleviated by the use of burst-allocation, as proposed by Ion Gaztañaga.

Too powerful to evolve

2008-09-18T07:11:00.000+02:00

Paragraph 1.4.8 of the C++ standard (draft copy) states that

A conforming implementation may have extensions (including additional library functions), provided they do not alter the behavior of any well-formed program.

It turns out that, strictly speaking, this unexpectedly bans many possible extensions an implementation of the C++ standard library might provide; C++ is such a powerful beast that a well-formed program might be made to depend on the non-existence of certain entities.

Function and function templates. The following program

‎#include <algorithm>
‎
‎using namespace std;
‎
‎int any_of(int count,int n1,...){/*...*/}
‎
‎int main()
‎{
‎  int n=any_of(2,1023,233);
‎}
‎

is valid C++03, but will fail with the upcoming C++0x revision of the standard, as this includes a function template any_of within <algorithm>:

‎template <class InputIterator, class Predicate>
‎bool any_of(InputIterator first, InputIterator last, Predicate pred);

that takes precedence over the user defined any_of due to the using namespace std directive.

Classes and class templates. Similarly, this C++03 valid program

‎#include <functional>
‎
‎using namespace std;
‎
‎template<typename T>
‎struct hash
‎{
‎  size_t operator()(const T& t){/*...*/}
‎};
‎
‎int main()
‎{
‎  hash<int>   h;
‎  size_t s=h(983265);
‎}
‎

will fail with C++0x due to the presence of an homonym hash in <functional>.

Member functions. The previous situations can be avoided if non-standard extensions are declared in separate headers or belong to a namespace other than std. These measures cannot be taken, however, when we are adding members to a standard class. The following is a totally artificial yet valid C++03 program:

‎#include <vector>
‎
‎template<class Vector>
‎int foo(char x){/*...*/}
‎
‎template <typename T,void (T::*)()> struct mfp{};
‎
‎template<class Vector>
‎void* foo(int x,mfp<Vector,&Vector::shrink_to_fit>* =0){/*...*/}
‎
‎int main()
‎{
‎  int x=foo<std::vector<int> >(0);
‎}

The so-called SFINAE rule prevents the second overload of foo to be selected because vector has no member function called shrink_to_fit; but it does in C++0x, rendering the program invalid (the second foo return type is not convertible to int).

Of course this is not a serious objection to the C++ standard text: 1.4.8 could always be rephrased to restrict it to reasonable well-formed programs, for some definition of "reasonable". But it can be instructive to reflect on how C++ has become so powerful (some would say complex) that such devious tricks can emerge (and, in some cases, be benefited from).

A causality calculus

2008-09-14T20:28:00.001+02:00

Let us design a first-order calculus including causality as non-truth-functional connective for which causality semantics is defined at the same level as the rest of the (truth-functional) connectives, according to the philosophical position stated in a previous entry. We assume we are given an infinite set of variables V and a set R of predicate symbols with associated arities.

Propositional constants: F (falsehood).

Connectives: → (only if), > (causal implication).

Quantifiers: V (for all).

Axioms and rules of inference: the usual ones for standard first-order calculus plus the following axiom schemas dealing with causality:

(C₁) p > q, q > r → p > r

(C₂) (p > p) → F

(C₃) p > q → (p and q)

(C₄) p > r → (p or q) > r

(C₅) p > (q and r) → p > q

These axioms deserve some comment: C₁ states the usual transitivity of causation, while C₂ says that no event causes itself; C₃ allows us to speak about causation only among events actually taking place, which makes this a "realist" semantics. C₄ and C₅ translate some usual properties of material implication to causal implication.

Semantics: a valuation v specifies several mappings:

From V to some set U.
From each n-ary predicate in R to some n-ary relation in U.
From each formula to the set {F,T}.

and obeys the following rules:

v(F) = F.
v(φ → ψ) = T iff v(φ) = F or v(ψ) = T.
v(R(t₁,...,t_n)) = T iff (v(t₁),...,v(t_n)) is in v(R).
v(Vx φ) = T iff v'(φ) = T for every valuation v' identical to v except possibly at x.
The relation on formulas of the calculus defined by v(φ > ψ) = T is consistent with axiom schemas C₁,..., C₅. In particular, this implies that the relation thus defined is a strict partial order.

As usual, a statement is valid if it is true for any valuation. It is easy to see that there are valuations with the properties described above: just take any valuation for standard first-order predicate calculus and augment it with v(φ > ψ) = F for all φ, ψ (we may call this a Humean valuation).

Introducing stable_vector

2008-09-10T04:08:00.000+02:00

We present stable_vector, a fully STL-compliant stable container that provides most of the features of std::vector except element contiguity. The underlying data structure has been discussed in a prior entry.

stable_vector.hpp. The code uses Boost.
test_stable_vector.cpp: Test suite for stable_vector.

In the following we are using the terminology of the C++98 standard.

General properties. stable_vector satisfies all the requirements of a container, a reversible container and a sequence and provides all the optional operations present in std::vector. Like std::vector, iterators are random access. stable_vector does not provide element contiguity; in exchange for this absence, the container is stable, i.e. references and iterators to an element of a stable_vector remain valid as long as the element is not erased, and an iterator that has been assigned the return value of end() always remain valid until the destruction of the associated stable_vector.

Operation complexity. The big-O complexities of stable_vector operations match exactly those of std::vector. In general, insertion/deletion is constant time at the end of the sequence and linear elsewhere. Unlike std::vector, stable_vector does not internally perform any value_type destruction, copy or assignment operations other than those exactly corresponding to the insertion of new elements or deletion of stored elements, which can sometimes compensate in terms of performance for the extra burden of doing more pointer manipulation and an additional allocation per element.

Exception safety. (according to Abrahams' terminology) As stable_vector does not internally copy elements around, some operations provide stronger exception safety guarantees than in std::vector:

operation	exception safety for `std::vector<T>`	exception safety for `stable_vector<T>`
`insert`	strong unless copy or assignment of `T` throw (basic)	strong
`erase`	no-throw unless copy or assignment of `T` throw (basic)	no-throw

Memory overhead. The C++ standard does not specifiy requirements on memory consumption, but virtually any implementation of std::vector has the same behavior wih respect to memory usage: the memory allocated by a vector v with n elements of type T is

m_v = c·e,

where c is v.capacity() and e is sizeof(T). c can be as low as n if the user has explicitly reserved the exact capacity required; otherwise, the average value c for a growing vector oscillates between 1.25·n and 1.5·n for typical resizing policies. For stable_vector, the memory usage is

m_sv = (c + 1)p + (n + 1)(e + p),

where p is the size of a pointer. We have c + 1 and n + 1 rather than c and n because a dummy node is needed at the end of the sequence. If we call f the capacity to size ratio c/n and assume that n is large enough, we have that

m_sv/m_v ≈ (fp + e + p)/fe.

So, stable_vector uses less memory than std::vector only when e > p and the capacity to size ratio exceeds a given threshold:

m_sv < m_v ↔ f > (e + p)/(e − p). (provided e > p)

This threshold approaches typical values of f below 1.5 when e > 5p; in a 32-bit architecture, when e > 20 bytes.

Summary. stable_vector is a drop-in replacement for std::vector providing stability of references and iterators, in exchange for missing element contiguity and also some performance and memory overhead. When the element objects are expensive to move around, the performance overhead can turn into a net performance gain for stable_vector if many middle insertions or deletions are performed or if resizing is very frequent. Similarly, if the elements are large there are situations when the memory used by stable_vector can actually be less than required by std::vector.