# How do I create a matrix of functions?

I'm trying to create a matrix in which the elements are functions of two variables, but I can't figure out how to do it. Is it possible? The only way I could figure out is to define the matrix in terms of two variables, then replace those variables with values later on. Here's my code:

row1 = {};
row2 = {};
For[ii = 1, ii <= Length[modes], ii++,
{l, m} = modes[[ii]];
G1lm = sglm[-2, l, m, \[Theta]] + (-1)^l sglm[-2, l, -m, \[Theta]];
G2lm = sglm[-2, l, m, \[Theta]] - (-1)^l sglm[-2, l, -m, \[Theta]];
AppendTo[row1, G1lm Cos[m \[Phi]]];
AppendTo[row1, G1lm Sin[m \[Phi]]];
AppendTo[row2, G2lm Sin[m \[Phi]]];
AppendTo[row2, -G2lm Cos[m \[Phi]]];
];
proj = {row1, row2};


Theta and Phi are my variables. Later on, when I want proj for a specific Theta and Phi I use replace:

pN = proj /. {\[Theta] -> a, \[Phi] -> b};


where a and b are some numbers. The problem with this solution is that if I use Theta and Phi for something else later on, say I set Theta = 0.5, then proj will permanently have Theta = 0.5, regardless of the replace I use. This happened to me: I was calling proj in a loop, but I kept getting the same results. After ripping my hair out for a few hours, I realized I had Theta set to a specific value in a block of code further up in my notebook. The simple solution is to simply not use Theta and Phi anywhere else in the notebook, but that doesn't seem very stable.

The root problem here is I'm using global variables to do something that should really be local to proj, but I don't know how to fix that. The ideal would be if I could somehow make proj a function of Theta and Phi so that later on I can simply do pN = proj[a,b]. I know how to make the Gnlm's functions of theta, but then I don't know how to make row1, row2, and by extension proj, functions of Theta (and Phi). Any help would be greatly appreciated.

-
What are modes and sglm? Otherwise I can only post a general answer on how to create a list of functions. Or rather, as it seems that is what you need, a function in two variables that creates a list. –  einbandi Mar 6 '13 at 17:11
modes is a 2x5 array of numbers --- specifically {{2,2},{2,1},{3,3},{4,4},{5,5}} --- and the sglm are the theta-dependent parts of the spin-weighted spherical harmonics. Honestly, the exact functional form of the sglm doesn't matter --- you could just use Sin[l theta] and Cos[l theta] if you wish. The point is each element of the matrix is a different function of the same two variables, theta and phi. –  ccap83 Mar 6 '13 at 17:20

You probably need something like this:

proj[modes_][θ_, ϕ_] :=
Table[
(sglm[-2, #1, #2, θ] +
sign (-1)^#2 sglm[-2, #1, -#2, θ])
If[sign == -1, {Cos[#2 ϕ],
Sin[#2 ϕ]}, {Sin[#2 ϕ], -Cos[#2 ϕ]}]&,
Transpose@modes],
{sign, {-1, 1}}
]


The list of modes is iterated through automatically, without the use of For, by applying MapThread to Transpose@modes.

The two lists you combine in proj differ only in the value of the variable sign which is implemented inside Table that goes through the values 1 and -1.

The function takes modes as its initial argument, and the angles are the other two arguments. You don't need to separate them like I did, that's just a matter of taste. The main thing here is how you use the underscore to declare the function variables. And the function itself, called proj, then returns the list.

In conclusion: functions can return lists, and For loops can be replaced by Map, MapThread or Table (etc.).

Edit

If you want to really pre-compute everything as far as possible given a fixed list modes, you can do the following:

Clear[proj, θ, ϕ]
modes = {{1, 0}, {1, 1}, {2, 2}};

proj[θ_, ϕ_] =
Table[Flatten[
sign (-1)^#1 sglm[-2, #1, -#2, θ]) If[
sign == 1, {Cos[#2 ϕ],
Sin[#2 ϕ]}, {Sin[#2 ϕ], -Cos[#2 ϕ]}] &,
Transpose@modes]], {sign, {1, -1}}]


(* ==> {{0,0,Cos[ϕ] (-sglm[-2,1,-1,θ]+sglm[-2,1,1,θ]),(-sglm[-2,1,-1,θ]+sglm[-2,1,1,θ]) Sin[ϕ],Cos[2 ϕ] (sglm[-2,2,-2,θ]+sglm[-2,2,2,θ]),(sglm[-2,2,-2,θ]+sglm[-2,2,2,θ]) Sin[2 ϕ]},{0,-2 sglm[-2,1,0,θ],(sglm[-2,1,-1,θ]+sglm[-2,1,1,θ]) Sin[ϕ],-Cos[ϕ] (sglm[-2,1,-1,θ]+sglm[-2,1,1,θ]),(-sglm[-2,2,-2,θ]+sglm[-2,2,2,θ]) Sin[2 ϕ],-Cos[2 ϕ] (-sglm[-2,2,-2,θ]+sglm[-2,2,2,θ])}} *)

This uses Set instead of SetDelayed to evaluate the modes and (if possible) sglm at the time you define the function proj. Then the MapThread won't have to be repeated when you call the function. It won't save significant time to do this unless parts of the sglm can in fact be evaluated for symbolic arguments.

-
Thanks! That seems to do the trick. It looks like this is building the table on the fly each time. Do you know if there is much computational cost in that? I'm wondering if it is any faster to pre-construct proj for the given modes, ala proj[\[Theta]_, \[Phi]_]:= proj[modes][\[Theta],\[Phi]]. I ask because proj is going inside FindMaximum. In that the number of modes always stays the same, just the angles change. –  ccap83 Mar 6 '13 at 17:32
Since the variables $\theta$ and $\phi$ are inside the sglmn, you would indeed have to re-do the sum for every new $\theta$ and $\phi$. That's why I don't see any significant room for improvement in terms of efficiency. But certainly you can leave out the modes argument if that stays fixed as a global variable. Then you don't even need to include it in my definition above, just make sure modes is defined once you call the function. –  Jens Mar 6 '13 at 17:42
As far as I see it, that would still mean the whole MapThread procedure is done every time the function is called. I think the OP wants to precompute everything, to get a function that really only inserts the arguments into the precomputed list. –  einbandi Mar 6 '13 at 17:46
@einbandi But that essentially costs nothing. –  Jens Mar 6 '13 at 18:01
Good point. I think it's interesting nevertheless, how one can compute complicated things and turn them to a function afterwards. Is there a standard way of replacing certain variables in an expression by Slots in order to create a function? –  einbandi Mar 6 '13 at 18:05

To build the table only when given (θ, ϕ) parameters I might use:

f1[θ_, ϕ_][{l_, m_}] :=
With[{
A = sglm[-2, l,  m, θ],
B = sglm[-2, l, -m, θ],
s = Sin[m ϕ],
c = Cos[m ϕ]
},
{{c, s}, {s, -c}} (A + (-1)^l {1, -1} B)
]

f2[θ_, ϕ_, modes : {{_,_} ..}] := Flatten[f1[θ, ϕ] /@ modes, {{2}, {1, 3}}]


Then for example:

f2[2.5, 7.2, {{2, 2}, {2, 1}, {3, 3}, {4, 4}, {5, 5}}]


In this construct all symbols are localized except for the function sglm and f1 and f2 themselves.
I could have put f1 in a Module inside f2 but I felt this presentation was cleaner.

-
Thanks, this works, but I'm finding it to be about 6x's slower than the replacement method I was using originally: Mean[Table[Timing[f2[0.5 Pi, 0.25 Pi, modes]][[1]], {i, 1, 1000}]] gives ~6e-4s, Mean[Table[ Timing[proj /. {\[Theta] -> 0.5 Pi, \[Phi] -> 0.25 Pi}][[1]], {i, 1, 1000}]] gives ~1e-4s. –  ccap83 Mar 6 '13 at 22:55
@ccap83 You did not seem to stress performance so I did not focus on that. If you provide the definition of your sglm function I can probably make my code faster, and perhaps yours too. –  Mr.Wizard Mar 9 '13 at 0:38

Using the question posted by @einbandi and answer posted by @acl here the general solution seems to be:

Module[{vars},
mat = {(*Some code that constructs a matrix out of expressions of vars*)}
fmat = Function @@{{vars}, mat};
]


Applying to my specific case, the following works:

Module[{\[Theta], \[Phi]},
row1 = {};
row2 = {};
For[ii = 1, ii <= Length[modes], ii++,
{l, m} = modes[[ii]];
G1lm = sglm[-2, l, m, \[Theta]] + (-1)^l sglm[-2, l, -m, \[Theta]];
G2lm = sglm[-2, l, m, \[Theta]] - (-1)^l sglm[-2, l, -m, \[Theta]];
AppendTo[row1, G1lm Cos[m \[Phi]]];
AppendTo[row1, G1lm Sin[m \[Phi]]];
AppendTo[row2, G2lm Sin[m \[Phi]]];
AppendTo[row2, -G2lm Cos[m \[Phi]]];
];
proj = Function @@ {{\[Theta], \[Phi]}, {row1, row2}};
]


Enclosing the construction of the matrix in a Module ensures that \[Theta] and \[Phi] are protected, and the last line turns proj into a function of the two variables. This method is also fast: it's about twice as fast as the replacement method I was originally using.

I'm sure there's much more concise ways to construct the matrix than the For loop I'm using here, but since this part of the code is only called once, any computational time savings are negligible.

-