2
$\begingroup$

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?

$\endgroup$
  • 2
    $\begingroup$ Check this question. I can't find Normalize and NormalDistribution $\endgroup$ – Dr. belisarius Feb 26 '16 at 7:48
  • 1
    $\begingroup$ Personally I also suggest you to read this answer: mathematica.stackexchange.com/a/104031/1871 $\endgroup$ – xzczd Feb 26 '16 at 15:37
  • $\begingroup$ @Dr.belisarius The Normalize and NormalDistribution are not the problem since it compiles if I remove MapThread $\endgroup$ – syhpphys Feb 26 '16 at 17:19
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Dr. belisarius Feb 26 '16 at 18:45
2
$\begingroup$

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]]
$\endgroup$
  • 1
    $\begingroup$ 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}. $\endgroup$ – syhpphys Feb 29 '16 at 22:35
  • $\begingroup$ 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? $\endgroup$ – syhpphys Feb 29 '16 at 22:46
  • 1
    $\begingroup$ @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. $\endgroup$ – xzczd Mar 1 '16 at 2:08
  • $\begingroup$ @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. $\endgroup$ – xzczd Mar 1 '16 at 2:10
  • 1
    $\begingroup$ @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 ^_^ $\endgroup$ – Alexey Golyshev Mar 1 '16 at 9:47
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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… $\endgroup$ – xzczd Feb 26 '16 at 15:34
  • 1
    $\begingroup$ "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. $\endgroup$ – xzczd Feb 26 '16 at 15:45
  • $\begingroup$ @xzczd Executed f = f[x] 20 times - no suppression. $RecursionLimit::reclim2: Recursion depth of 1024 exceeded during evaluation of f[x]. $\endgroup$ – Alexey Golyshev Feb 26 '16 at 15:52
  • $\begingroup$ 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.stack.imgur.com/tOlHH.png $\endgroup$ – xzczd Feb 26 '16 at 16:01
  • $\begingroup$ @xzczd v10.3.1 fastswf.com/jgyvaTM $\endgroup$ – Alexey Golyshev Feb 26 '16 at 16:06
1
$\begingroup$

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

Mathematica graphics

$\endgroup$
  • $\begingroup$ It works, but it is slower compared to the original version when not compiled; speed is the only reason I am compiling. $\endgroup$ – syhpphys Feb 26 '16 at 22:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.