Friday, November 28, 2008

More on price gas hysteresis

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 + mx = 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 + mx = 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.

Monday, November 24, 2008

Optimizing red-black tree color bits

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_;

‎ T& value(){return value_;}
‎ const T& value()const{return value_;}

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

  1. This is requested by the user via ColorBitCompression.
  2. There is an unsigned integer with which pointers can be represented.
  3. 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_;

‎ 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_;

‎ 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.

Thursday, November 20, 2008

Infinitely many voters

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={v1,...,vN}: we "extend" them by assigning to each vi an enumerable set of fictitious voters vi1, vi2,... 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 vi is obtained and assigned to the associated fictitious voters vij; 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 : VV' is not the entire powerset of V', but only a finite subset of it (with cardinality 2N); 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 vi, who is the resulting dictator for the so devised voting system.

Sunday, November 16, 2008

State semantics

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 BT whose internal state is represented by values of T. We replace T::f with the following immutable method:
‎BT BT::f(args)
‎ T state=this->state;
‎ state.f(args);
‎ return BT(state);
so that the expression
‎T x;
can be rewritten with BT as
‎BT x;
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 ST whose internal state is represented by values of T, and implement the associated method f as follows:
‎void ST::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.

Wednesday, November 12, 2008

Gas price hysteresis

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 + mx = 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 + mx = 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.

Saturday, November 8, 2008

Multiattribute querying with Boost.MultiIndex

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 columni = valuei 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

‎// x1==5 && x4==40 && x2==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 -->
‎x1==v1 && x3==v3 && x4==v4 -->
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:
‎ 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, \
‎{ \
‎ 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< \
‎ > fields_type; \
‎ \
‎ const index_type& i=c.get<tag>(); \
‎ boost::tuple< \
‎ > t=boost::make_tuple( \
‎ ); \
‎ 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; \

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.

Tuesday, November 4, 2008

Filling up at dawn

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

(mcold mhot)/mhot = 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) = ΔTsurface·ez/D,
D = √(2K/ω),

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

1.9·10-12ΔTsurface ≤ ΔTstation tank ≤ 1.4·10-6ΔTsurface,

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

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

we have

7.3 °C ≤ ΔTstation tank ≤ 14.8 °C,

which translates to

0.7% ≤ (mcold mhot)/mhot t ≤ 1.4%.

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

Friday, October 31, 2008

C preprocessor tail recursion

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) \
‎ n, \
‎ 1)
The expression FACTORIAL(5) does not expand to 120, however, but to

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

The last line of the snippet above does not enter into infinite recursion, but simply causes the expansion rules ABA 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 ff' 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:
‎ 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) \

‎#define RECURSE_PRED(d,state) \

‎#define RECURSE_OP(d,state) \

‎#define RECURSE_BODY(state) \

‎#define RECURSE_RES(state) \

‎#define STOP (0)

‎#define ITERATE (1)
which allows us to implement our factorial macro in a reasonably straightforward manner:

‎#define FACTORIAL_BODY(args) \

‎#define FACTORIAL_BODY_IMPL(n,res) \
‎ n, \
‎ 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.

Monday, October 27, 2008

Law-abiding miracles

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.

Thursday, October 23, 2008

Circles, ellipses, assistants and employees

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) = [t0,t1] 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 R2: 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.

Monday, October 20, 2008

Multicolumn querying

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 (c1,...,cn), sorting is performed according to the lexicographical order associated to the n columns, with ci taking precedence over cj 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 (c1,...,cn) serves to accelerate queries of the form

c1 = v1 AND c2 = v2 AND ··· AND ci = vi

for all in. 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 2n − 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(2n/√n). The mathematical formulation of this problem has been studied in a previous entry. Also, an algorithm for generating the index set is provided.

Thursday, October 16, 2008

On necessity in McGinn

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.

Sunday, October 12, 2008

Points inside a polygon and the Residue Theorem

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 Resz=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 = (p1,p2), (p2,p3), ... , (pN−1,pN), (pN,p1),

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

P (1/z) dz = (pi,pi+1) (1/z) dz =
= (ln pi+1 − ln pi).

But now, this sum seems to be unconditionally null:

(ln pi+1 − ln pi) =
= ln p2 − ln p1 + ln p3 − ln p2 + ln p4 − ln p3 + ... + ln p1 − ln pN = 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

(pi,pi+1) (1/z) dz = ln pi+1 − ln pi

only holds if the derivative of ln z is 1/z along all the segment (pi,pi+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 (pi,pi+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 (pi,pi+1) does not intersect the new cut. It is not hard to see that doing so yields

(pi,pi+1) (1/z) dz = ln pi+1 − ln pi + σ2πi,

where ln denotes again the principal branch of the complex logarithm and σ is +1 or −1 depending on whether (pi,pi+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 =(sn)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.

Wednesday, October 8, 2008

Generating permutation covers: part II

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 S0 U ... U SM, where Si 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 S0 U ... U Sm, 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 S0 U ... U Sm U Sm + 1. We state without proof the following

Lemma. If mN/2, a miminal tuple cover of S0 U ... U Sm 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 Sm + 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 Sm must be covered by some tuple of T. So the extension mapping T T' relates elements of Sm to elements of Sm + 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 Sm is the source of (Nm) arrows, whereas m + 1 arrows arrive at each element of Sm + 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 Sm + 1 an arbitrary arrow arriving at it there is the possibility that some elements of Sm are not visited by any of the selected arrows. The following is a balanced selection criterion: given a fixed injection f : SmSm + 1 and function g : Sm + 1Sm, both compatible with the inclusion relationship between Sm and Sm + 1 (that is, X is a subset of f(X) and Y a superset of g(Y)), we select an arrow XY iff

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

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 Sm + 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 Sm to Sm + 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 χ : Sm + 1 → {true,false} that checks whether a given Y belongs to f(Sm). We adopt the following definition for g:

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

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

Τ' ← Ø
for every τ in Τ
····X ← range(τ)
····af(X) − X {the one element added to X by f}
····Τ'Τ' U {τ·a}
····for i = max(X) + 1 to N − 1
········YX 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 XY arrows where X = g(Y), which is maximally efficient and saves us the need to explicity compute g. The complete Τ0 = Ø → Τ1 → ··· → Τ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 {τ}
····X ← range(τ)
····af(X) − X
····tuple-cover(m + 1,N,τ·a,Τ)
····for i = max(X) + 1 to N − 1
········YX 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 : SMSM such that Xd(X) = Ø for all X in SM, 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}.

if N is even then
····N'N − 1
end if
Σ ← Ø
Τ ← Ø
for every τ in Τ
····find τ' in Τ with range(τ') = d(range(τ))
····aS − 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
2(ab), (ba)
3(acb), (bac), (cba)

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