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$
3
  • 2
    $\begingroup$ Check this question. I can't find Normalize and NormalDistribution $\endgroup$ Commented Feb 26, 2016 at 7:48
  • 1
    $\begingroup$ Personally I also suggest you to read this answer: mathematica.stackexchange.com/a/104031/1871 $\endgroup$
    – xzczd
    Commented Feb 26, 2016 at 15:37
  • $\begingroup$ @Dr.belisarius The Normalize and NormalDistribution are not the problem since it compiles if I remove MapThread $\endgroup$
    – syhpphys
    Commented Feb 26, 2016 at 17:19

3 Answers 3

3
$\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$
6
  • 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
    Commented Feb 29, 2016 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
    Commented Feb 29, 2016 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
    Commented Mar 1, 2016 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
    Commented Mar 1, 2016 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$ Commented Mar 1, 2016 at 9:47
3
$\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$
6
  • $\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
    Commented Feb 26, 2016 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
    Commented Feb 26, 2016 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$ Commented Feb 26, 2016 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.sstatic.net/tOlHH.png $\endgroup$
    – xzczd
    Commented Feb 26, 2016 at 16:01
  • $\begingroup$ @xzczd v10.3.1 fastswf.com/jgyvaTM $\endgroup$ Commented Feb 26, 2016 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$
1
  • $\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
    Commented Feb 26, 2016 at 22:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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