type inference changing dimensionality

I use a technique called ‘type reconstruction’, which is a combination of function specialization and type inference to find the best matching type in this case. I will have to look in more detail to your example why the resulting type isn’t mat.
In any case, although I would recommend to avoid it in general, type casting can still be performed using the cast function. This function is mostly there for correctly handling of generic conversions in a function with generic parameters.
For example:

B = cast(0.0, T)  %where T=scalar or T=cscalar
y = cast(bar, mat)

OK, I have checked. The generic narrowing from higher dimensional to lower dimensional structures (e.g. cube to mat) has not been implemented yet in this context. Here you provide a very interesting use case. Currently, the compiler will only apply the type reconstruction when at least one of the output types of the function is unspecified (so this is the generic case you mention yourself).
There is a trick to test how the type reconstruction sub-system is working: just pull in the following line:

 print $typerecon(generic_function, type(A,mat))
 % Result is a function type: [cube->cube] or [mat->mat]

The second argument for $typerecon can be a logic expression (e.g. type(A,mat) && B==4 && C==8 etc.)
The interesting part is that the type narrowing boils down to an extension of the evidence handling functions in the truth value system. Which means that it may have a (positive) effect on other internal components as well. I will check how difficult it is to make this extension.

Good news: it works! The extension only required minor changes and did not impact the regression tests, apart from one small issue that I need to look at later (generic iteration over multidimensional structures).
Speaking about other effects: this technique can be used to reduce the dimension of all loops within a function (i.e., loop flattening but then by specialization). In the following example:

function y = binary_op[T](a : cube[T], b : cube[T], op : [__device__ (T,T)->T]) 
    assert(size(a)==size(b))
    y = copy(a)
    function [] = __kernel__ kernel (a, b, y, op, pos)
        y[pos] = op(a[pos], b[pos])
    end

    parallel_do(size(a),a,b,y,op,kernel)
end

The specialization:

binary_op_vec = $specialize(binary_op, type(a, vec) && type(b, vec))

will then generate code in which the length of the pos argument in the inner function is reduced (1 instead of 3):

% Specialization of function 'binary_op' with constraints '(type(a,vec)&&type(b,vec))'
function y=binary_op_vec(a:vec,b:vec, op:[__device__ (scalar,scalar)->scalar])
    assert((size(a)==size(b)),"Note: passed compiler check.")
    y=copy(a)

    % Specialization of function 'kernel' with constraints '(((type(a,vec)&&type(b,vec))&&type(y,vec))&&type(pos,int))'
    function [] = __kernel__ kernel(a:vec'const,b:vec'const,y:vec,op:[__device__ (scalar,scalar)->scalar],pos:int)
       y[pos] = op(a[pos], b[pos])       
    end
    parallel_do(size(a),a,b,y,op,kernel)
end

Of course, when the indexing is not as trivial as here, the function specialization may fail resulting in a compiler error (e.g. indexing with too many indices) – but this would then be as expected. This does not affect the type reconstruction, so the above extension is in any case a good thing.

Thanks for the fix!
The other effects seem exciting too, the loop flattening in particular. But is the indexing breakage intended? Why shouldn’t indexing a specialized lower-dimensional type work?
In Julia, indexing with too many indices never fails (if they are in-bounds, of course):

function test(x::Array)
    println(x[1, 1])
end

test(randn(2, 2))
test(randn(2))

This allows the compiler to specialize test for a lower-dimensional input without risking to break indexing.
In Quasar, this doesn’t seem to work (or at least not reliably):

a = randn(2, 2)
print type(a)           % mat
print size(a, 0..10)    % [2, 2, 1...]
print a[0, 0]           % works as expected
print a[0, 0, 0]        % works?
print a[0, 0, 0, 0]     % out of range?

This is a good suggestion and could potentially solve the dimensionality reduction problem with indices that I mentioned above. The only thing is that it moves the problem to the run-time which may lead to out of bounds errors there (these “bugs” are generally more difficult to solve for the user). Therefore, it would be good that the compiler checks the bounds as well (i.e. for all singleton dimensions the index should be zero), resulting in an appropriate warning/error.
I also need to check the indexing cases that you mentioned (a[0,0,0], a[0,0,0,0]) and ask around if everyone would like this Julia indexing feature (could have unintended side effects when one forgets to remove an index when converting from cube to mat; this is also why it is not there yet).