On Wednesday, May 20, 2015 at 2:22:59 PM UTC-4, Júlio Hoffimann wrote:
>
> David, thank you very much for the detailed explanation, will read through
> it.
>
> Steven, sorry for not explaining the problem clearly. My issue is
> basically trying to figure out a way to loop over all possible combinations
> of exponents (e1, e2, e3) for which the sum is less or equal to N and
> generate a vector with the terms Xj[1]^e1*Xj[2]^e2*Xj[3]^e3.
>
> The j-th column of the resulting matrix contains the terms of this
> polynomial evaluated at the j-th column of X.
>
Thanks for the explanation.
If N is fixed (or you have a small number of N, known in advance), and you
want to evaluate this for lots of different 3xM matrices, and performance
is really critical, it is not crazy to use metaprogramming to generate an
optimized inlined version of this.
Here an example that I think does what you want:
macro generate_vanderish(N)
e123 = Array(Tuple{Int,Int,Int,Int},0) # (index,e1,e2,e3) exponents
i = 0
for e1=0:N
for e2=0:N-e1
i += 1
push!(e123,(i,e1,e2,N-e1-e2))
end
end
M = length(e123)
vanderish! = esc(symbol(string("vanderish",N,"!")))
vanderish = esc(symbol(string("vanderish",N,)))
quote
function $vanderish!(Y, X)
size(X,1) != 3 && error("X must be 3xM matrix")
size(Y,2) != size(X,2) && error("mismatched #cols in Y and X")
size(Y,1) != $M && error("Y must have ", $M, " rows")
for j = 1:size(X,2)
x1 = X[1,j]
x2 = X[2,j]
x3 = X[3,j]
$(Expr(:block,
[:(Y[$i,j] = x1^$e1 * x2^$e2 * x3^$e3)
for (i,e1,e2,e3) in e123]...))
end
return Y
end
function $vanderish(X)
Y = similar(X, $M, size(X,2))
$vanderish!(Y, X)
end
end
end
For example, try
@generate_vanderish(3)
which defines a function vanderish3(X) that operates on a 3xM matrix X and
outputs your "Vandermonde-like" matrix. It also defines an in-place
function vanderish3!(Y,X) in which you supply a pre-allocated output matrix
Y.
A good way to understand how this works is to change the word "macro" to
"function" above, and just call generate_vanderish(3) to see the expression
that it outputs.
If you have to handle matrices X with lots of different numbers of rows
(not just 3), then code generation makes less sense.
