3
$\begingroup$

This is an extension of my earlier questions Faster use of Condition for a large array (see the example) and How to use Compile to generate an $n\times n$ array using $n$ vectors . In this case each point is characterised by a pair of indices.

Consider this Array.

n1 = n2 = 30;
n = 9;
data = Flatten[Table[{RandomInteger[{1, n}, 2], {1. i, 1. j}},
{i, n1}, {j, n2}],1];

So each element of data is {{Integer,Integer},{Real,Real}}. Then I define 3 matrices.

P = RandomReal[{0, 1}, {n, n}];
Q = RandomReal[{0, 1}, {n, n}];
R = RandomReal[{0, 1}, {n, n}];

And I have to create another matrix from data like this

m1 = Table[
{s1, t1} = x[[1]]; {s2, t2} = y[[1]];
dr = x[[2]] - y[[2]];
Which[
 Abs[Norm[dr] - 0.] < 0.001, R[[s1, t1]],
 Abs[Norm[dr] - 1.] < 0.001, Q[[Abs[s1 - s2] + 1, Abs[t1 - t2] + 1]],
 Abs[Norm[dr] - Sqrt[2]] < 0.001, P[[Abs[s1 - s2] + 1, Abs[t1 - t2] + 1]],
 True, 0],
{x, data}, {y, data}]; // AbsoluteTiming

{17.569614, Null}

That is, depending on what is distance between $i^{th}$ and $j^{th}$ point, the $ij^{th}$ element of m1 will be given by an element of P,Q or R matrix.

As much I know, in such cases Compile gives faster result. So I try in this way (based on the answer to a similar problem by Pickett)

TensorRank[data]

3

m = Compile[{{sys, _Real, 3}},
Table[
Assuming[{s1, t1, s2, t2} \[Element] Integers,
{{s1, t1} = IntegerPart@x[[1]], {s2, t2} = IntegerPart@y[[1]]}];
dr = x[[2]] - y[[2]];
Which[
 Abs[Norm[dr] - 0.] < 0.001, R[[s1, t1]],
 Abs[Norm[dr] - 1.] < 0.001, Q[[Abs[s1 - s2] + 1, Abs[t1 - t2] + 1]],
 Abs[Norm[dr] - Sqrt[2]] < 0.001, P[[Abs[s1 - s2] + 1, Abs[t1 - t2] + 1]],
 True, 0],
{x, sys}, {y, sys}],
CompilationTarget -> "C"];

m2 = m[data]; // AbsoluteTiming

{48.353723, Null}

It is almost same time with "RuntimeOptions" -> "Speed".

However if I remove Assuming, so it treats s1,t1,s2,t2 as Reals, it produces an error message

CompiledFunction::cfse: Compiled expression Compile`GetElement[{7,5},1] should be a machine-size integer. >>

CompiledFunction::cfexe: Could not complete external evaluation; proceeding with uncompiled evaluation. >>

and then it takes same time as the Table (~18sec).

What would be the correct and faster way to do this?

$\endgroup$
3
  • $\begingroup$ Try this: Needs["CompiledFunctionTools`"]; CompilePrint[m]. You will see many calls to MainEvaluate in your generated code that indicate calls to the main Mathematica evaluator for functions / expressions that cannot be compiled. You will want to take a look at How to Compile effectively? as well. $\endgroup$
    – MarcoB
    Aug 13, 2015 at 14:29
  • $\begingroup$ Could you describe in words what you want to accomplish? I can't quite follow your implementation. $\endgroup$
    – MarcoB
    Aug 13, 2015 at 14:48
  • $\begingroup$ Thanks @MarcoB. I know this question and the suggestions, they are indeed quite helpful. Regarding what I am trying to do - I have a set of n points characterized by two indices (first two integers) and coordinates (second pair of real numbers). Now I am trying to construct an nxn matrix whose ij-th element is determined by the coordinate and index of i-th and j-th point. Depending on the distance between ith and jth point (dr= 0.,1.,Sqrt[2.]) the ij-th element of m matrix will be given by an element of P,Q or R matrix, as described by the Which. Is it clear now? $\endgroup$
    – Sumit
    Aug 13, 2015 at 20:01

1 Answer 1

4
$\begingroup$

I'm not sure how to answer this question without sounding discouraging. I cannot go into every detail but I still want to help you. Basically, you did several things wrong:

  • using global variables inside Compile which makes that the compiled function calls back to the Mathematica Kernel to ask for instance what the definition of R is. Compiled functions need to be completely self-contained (until you are knowledgeable enough to know when to break this rule). This means, every definition, variable, etc that is not given in the arguments to Compile needs to be injected. Marco already gave you the correct links for further reading.

  • local variables need to put inside a Module (or With if appropriate)

  • not all functions of Mathematica can be used inside a compiled construct. You can check the list of compilable functions. When a function like Integrate cannot be compiled, then (again) the compiled function calls back to the Kernel asking him to evaluate certain parts.

  • type inference, which means how the compiler decides whether dr for instance is a real or integer variable is especially difficult inside Compile. The easiest way is to write s1=0 in the Module declaration to indicate s1 is an integer. Or s1=0.0 to make it a real. This issue was discussed here in on a different occasion.

  • Use CompilePrint from the <<CompiledFunctionTools` package to check if there are MainEvaluate calls in your compiled code that indicate that the Kernel is called.

In conclusion one can say, this all takes a lot of experience and it is never (except for the most simple cases) a solution to just write Compile around an expression. That being said, I have corrected your code. You can get my complete notebook by evaluating :

Import["http://goo.gl/NaH6rM"]["https://i.stack.imgur.com/S2KvJ.png"]

or you copy the pieces when reading the explanations below:

n1 = n2 = 30;
n = 9;
data = Flatten[
   Table[{RandomInteger[{1, n}, 2], {1. i, 1. j}}, {i, n1}, {j, n2}], 
   1];

P = RandomReal[{0, 1}, {n, n}];
Q = RandomReal[{0, 1}, {n, n}];
R = RandomReal[{0, 1}, {n, n}];

With[{P = P, Q = Q, R = R},
 fc1 = Compile[{{sys, _Integer, 3}},
    Module[{s1 = 0, t1 = 0, s2 = 0, t2 = 0, dr},
     Table[
      {s1, t1} = x[[1]];
      {s2, t2} = y[[1]];
      dr = Sqrt[#.#] &[x[[2]] - y[[2]]];
      Which[
       Abs[dr - 0.] < 0.001,
       R[[s1, t1]],
       Abs[dr - 1.] < 0.001,
       Q[[Abs[s1 - s2] + 1, Abs[t1 - t2] + 1]],
       Abs[dr - Sqrt[2]] < 0.001,
       P[[Abs[s1 - s2] + 1, Abs[t1 - t2] + 1]],
       True, 0.0], {x, sys}, {y, sys}]
     ], CompilationTarget -> "C"];
 ]

This function gives the same result as your Table call but only takes

fc1[data]; // AbsoluteTiming
(* {0.218867, Null} *)

0.2 seconds to run. There is still room for improvement. Your compiled function consists basically of a Table having two iterators x and y. What you could do is re-use the code above and adapt it so that you create a function that

iterates only over y and does the iteration over x in parallel.

To understand how this works, please see my answer in Use Compile to improve performance

With[{P = P, Q = Q, R = R},
 f2 = Compile[{{sys, _Integer, 3}, {x, _Integer, 2}},
    Module[{s1 = 0, t1 = 0, s2 = 0, t2 = 0, dr},
     Table[
      {s1, t1} = x[[1]];
      {s2, t2} = y[[1]];
      dr = Sqrt[#.#] &[x[[2]] - y[[2]]];
      Which[
       Abs[dr - 0.] < 0.001,
       R[[s1, t1]],
       Abs[dr - 1.] < 0.001,
       Q[[Abs[s1 - s2] + 1, Abs[t1 - t2] + 1]],
       Abs[dr - Sqrt[2]] < 0.001,
       P[[Abs[s1 - s2] + 1, Abs[t1 - t2] + 1]],
       True, 0.0], {y, sys}]
     ],
    CompilationTarget -> "C",
    Parallelization -> True,
    RuntimeAttributes -> {Listable},
    RuntimeOptions -> "Speed"
    ];
 ]

you have to call this function with fc2[data, data] and if you really put some effort in understanding what happens, then you will see why and how it works.

This brings an additional factor of 2

<< GeneralUtilities`

AccurateTiming[fc1[data]]/AccurateTiming[fc2[data, data]]
(* 2.06425 *)
$\endgroup$
3
  • $\begingroup$ “This issue was [discussed here] in on a different occasion.” You missed a link :) $\endgroup$
    – xzczd
    Aug 18, 2015 at 2:35
  • 1
    $\begingroup$ @xzczd That was left as exercise for the attentive reader. Fixed. $\endgroup$
    – halirutan
    Aug 18, 2015 at 2:43
  • $\begingroup$ Thanks @halirutan. That is very helpful for a beginner like me. $\endgroup$
    – Sumit
    Aug 18, 2015 at 6:33

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.