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?

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ Check this question. I can't find Normalize and NormalDistribution $\endgroup$
    – Dr. belisarius
    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
    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
    Feb 26, 2016 at 17:19

3 Answers 3

Reset to default
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]]
Improve this answer
$\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
    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
    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
    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
    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$
    – Alexey Golyshev
    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.

Improve this answer
$\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
    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
    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$
    – Alexey Golyshev
    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.stack.imgur.com/tOlHH.png $\endgroup$
    – xzczd
    Feb 26, 2016 at 16:01
  • $\begingroup$ @xzczd v10.3.1 fastswf.com/jgyvaTM $\endgroup$
    – Alexey Golyshev
    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

Improve this answer
$\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
    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.