Q:

What does the c underscore expression `c_` do exactly?

It seems to be some kind of horizontal concatenation, but I could not find any documentation online. Here a minimal working example: 
In [1]: from numpy import c_
In [2]: a = ones(4)
In [3]: b = zeros((4,10))    
In [4]: c_[a,b]
Out[4]: 
array([[ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])

A:

Use IPython's ? syntax to get more information:
In [2]: c_?
Type:       CClass
Base Class: <class 'numpy.lib.index_tricks.CClass'>
String Form:<numpy.lib.index_tricks.CClass object at 0x9a848cc>
Namespace:  Interactive
Length:     0
File:       /usr/lib/python2.7/dist-packages/numpy/lib/index_tricks.py
Docstring:
Translates slice objects to concatenation along the second axis.

This is short-hand for ``np.r_['-1,2,0', index expression]``, which is
useful because of its common occurrence. In particular, arrays will be
stacked along their last axis after being upgraded to at least 2-D with
1's post-pended to the shape (column vectors made out of 1-D arrays).

For detailed documentation, see `r_`.

Examples
--------
>>> np.c_[np.array([[1,2,3]]), 0, 0, np.array([[4,5,6]])]
array([[1, 2, 3, 0, 0, 4, 5, 6]])

A:

It took me a lot of time to understand but it seems I finally got it.
All you have to do is add along second axis.
let's take :
np.c_[np.array([1,2,3]), np.array([4,5,6])]

But there isn't second axis. So we mentally add one.
so shape of both array  becomes (3,1).
So resultant shape would be (3,1+1) which is (3,2).
which is the shape of result - 
array([[1, 4],
       [2, 5],
       [3, 6]])

Another ex.:
np.c_[np.array([[1,2,3]]), 0, 0, np.array([[4,5,6]])]

shapes: 
np.array([[1,2,3]]) = 1,3
np.array([[4,5,6]]) = 1,3
0 so we can think of it as [[0]] = 1,1
So result 1,3+1+1+3 = 1,8
which is the shape of result : array([[1, 2, 3, 0, 0, 4, 5, 6]])

A:

I would explain this as follow.  It concats your first array into the last dimension (axis) of your last array in the function.
For example:
# both are 2 dimensional array
a = array([[1, 2, 3], [4, 5, 6]])
b = array([[7, 8, 9], [10, 11, 12]])

Now, let's take a look at np.c_(a, b):
First, let's look at the shape:
The shape of both a and b are (2, 3).  Concating a (2, 3) into the last axis of b (3), while keeping other axises unchanged (1) will become 
(2, 3 + 3) = (2, 6)

That's the new shape.
Now, let's look at the result:
In b, the 2 items in the last axis are:
1st: [7, 8, 9]
2nd: [10, 11, 12]

Adding a to it means:
1st item: [1,2,3] + [7,8,9] = [1,2,3,7,8,9]
2nd item: [4,5,6] + [10,11,12] = [4,5,6,10,11,12]

So, the result is
[
  [1,2,3,7,8,9],
  [4,5,6,10,11,12]
]

It's shape is (2, 6)

