I have the following code:

s[x_Real, y_Real, z_Real, w_Real] = {{x + I y, - z + I w}, { z + I w, x - I y}};
ssamp = Compile[{{n, _Integer}, {nd, _Integer}}, MapThread[s, Transpose[Map[
Normalize, RandomVariate[NormalDistribution[0, 1], {n, nd, 4}],{2}], {3, 2, 1}] , 2]]


The compile fails with error messages

• MapThread::list: List expected at position 2
• Compile::cprank: Compile cannot determine the rank of the result tensor.

If I remove the MapThread, the function compiles. Any suggestions for getting it to compile?

• Check this question. I can't find Normalize and NormalDistribution Commented Feb 26, 2016 at 7:48
• Personally I also suggest you to read this answer: mathematica.stackexchange.com/a/104031/1871 Commented Feb 26, 2016 at 15:37
• @Dr.belisarius The Normalize and NormalDistribution are not the problem since it compiles if I remove MapThread Commented Feb 26, 2016 at 17:19

Here's a fully compiled version, principles applied in modifying the code has been mentioned in the links above:

ssamp2 = Compile[{{n, _Integer}, {nd, _Integer}},
Partition[Function[xyzw, {{xyzw[[1]] + I xyzw[[2]],
-xyzw[[3]] + I xyzw[[4]]},
{xyzw[[3]] + I xyzw[[4]],
xyzw[[1]] - I xyzw[[2]]}}]@#/Sqrt@Total[#^2] & /@
RandomVariate[NormalDistribution[0, 1], {n nd, 4}], n]]

• Now that you've written it so nicely, it is clear that we can throw away the Flatten and Transpose and just replace the {n,nd,4} with {n*nd,4}. Commented Feb 29, 2016 at 22:35
• The key thing that makes this work seems to be to normalize the vectors first, then pass them to the function defined within the Function[ , ]. May I ask how you figured this out? Commented Feb 29, 2016 at 22:46
• @syhpphys I think the 2 key points here are 1. Normalize can not be compiled so one needs to replace it with a compiled equivalent i.e. #/Sqrt@Total[#^2] & . 2. function defined with pattern-matching cannot be compiled (see the 3rd rule here?), so one needs to replace it with a pure function, which can usually be compiled if it's formed by compilable function. Commented Mar 1, 2016 at 2:08
• @syhpphys "we can throw away the Flatten and Transpose and just replace the {n,nd,4} with {n*nd,4}" You are right, I forgot to take the code inside RandomVariate into account, thanks for pointing out, edited. Commented Mar 1, 2016 at 2:10
• @xzczd +1 to your answer. In comments you wrote: "function defined with pattern-matching cannot be compiled". Very good comment, not only this phrase. Also +1. But in the next Mathematica it will be possible: youtu.be/… A lot of cool things expected ^_^ Commented Mar 1, 2016 at 9:47

You can replace MapThread[func,list,2] with func/@Transpose[Flatten/@list].

s[{x_Real,y_Real,z_Real,w_Real}]:={{x+I y,-z+I w},{z+I w,x-I y}};

ssamp = Compile[{{n, _Integer}, {nd, _Integer}},
s /@
Transpose[
Flatten /@
Transpose[
Map[Normalize,
RandomVariate[NormalDistribution[0, 1], {n, nd, 4}], {2}], {3, 2, 1}]
]
]


When you evaluate this compiled function, first 3 times and only 3 times (it's very strange!) you will see such error message.

ssamp[2, 2]


CompiledFunction::cfex: Could not complete external evaluation at instruction 36; proceeding with uncompiled evaluation. >>

But CompilePrint looks okay.

Needs["CompiledFunctionTools"]
CompilePrint@ssamp

"
2 arguments
17 Integer registers
2 Real registers
7 Tensor registers
Underflow checking off
Overflow checking off
Integer overflow checking on
RuntimeAttributes -> {}

I0 = A1
I1 = A2
I3 = 0
I2 = 4
T(I1)6 = {3, 2, 1}
I14 = -1
I4 = 1
I16 = 3
Result = T(R2)0

1   T(I1)1 = {I0, I1, I2}
2   R0 = I3
3   R1 = I4
4   T(R3)3 = RandomNormal[ R0, R1, T(I1)1]]
5   I8 = Length[ T(R3)3]
6   T(R2)0 = Part[ T(R3)3, I4]
7   I11 = Length[ T(R2)0]
8   I15 = I14
9   T(R3)0 = Table[ I8, I11, I15]
10  I12 = I3
11  goto 18
12  I13 = I3
13  goto 17
14  T(R1)4 = GetElement[ T(R3)3, I12, I13]
15  T(R1)5 = MainEvaluate[ Hold[Normalize][ T(R1)4]]
16  Element[ T(R2)0, I15] = T(R1)5
17  if[ ++ I13 <= I11] goto 14
18  if[ ++ I12 <= I8] goto 12
19  T(R3)3 = Transpose[ T(R3)0, T(I1)6, I16]]
20  I6 = Length[ T(R3)3]
21  I9 = I14
22  T(R2)0 = Table[ I6, I9]
23  I7 = I3
24  goto 28
25  T(R2)5 = GetElement[ T(R3)3, I7]
26  T(R1)4 = Flatten[ T(R2)5, I4]]
27  Element[ T(R2)0, I9] = T(R1)4
28  if[ ++ I7 <= I6] goto 25
29  T(R2)3 = Transpose[ T(R2)0]]
30  I5 = Length[ T(R2)3]
31  I6 = I14
32  T(R2)0 = Table[ I5, I6]
33  I9 = I3
34  goto 38
35  T(R1)4 = GetElement[ T(R2)3, I9]
36  T(R1)5 = MainEvaluate[ Hold[s][ T(R1)4]]
37  Element[ T(R2)0, I6] = T(R1)5
38  if[ ++ I9 <= I5] goto 35
39  Return
"


I hope this helps.

• Well, I don't think the 15 T(R1)5 = MainEvaluate[ Hold[Normalize][ T(R1)4]] and 36 T(R1)5 = MainEvaluate[ Hold[s][ T(R1)4]]  lines look OK… Commented Feb 26, 2016 at 15:34
• "When you evaluate this compiled function, first 3 times and only 3 times (it's very strange!) you will see such error message" it's because the warning is suppressed then. (I understand it as a silent General::stop. ) To see a similar behavior, execute f = f[x] for 3 times. Commented Feb 26, 2016 at 15:45
• @xzczd Executed f = f[x] 20 times - no suppression. \$RecursionLimit::reclim2: Recursion depth of 1024 exceeded during evaluation of f[x]. Commented Feb 26, 2016 at 15:52
• What's your version? I observed this behavior in v8.0.4 and v9.0.1. Here's the screenshot for v8.0.4: i.sstatic.net/tOlHH.png Commented Feb 26, 2016 at 16:01
• @xzczd v10.3.1 fastswf.com/jgyvaTM Commented Feb 26, 2016 at 16:06

The following seems to work (at least on v9):

s[{x_Real, y_Real, z_Real, w_Real}] := {{x + I y, -z + I w}, {z + I w, x - I y}}

ssamp = Compile[{{n, _Integer}, {nd, _Integer}},
Module[{u1},
u1 = Transpose[
Map[Normalize, RandomVariate[NormalDistribution[0, 1], {n, nd, 4}], {2}],
{3, 2, 1}];
s /@ Transpose[Flatten /@ u1]],
{{s[__], _Complex, 2}},
CompilationOptions -> {"InlineExternalDefinitions" -> True}];

ssamp[2, 2] // Grid
`

• It works, but it is slower compared to the original version when not compiled; speed is the only reason I am compiling. Commented Feb 26, 2016 at 22:03