6
$\begingroup$

I am trying to calculate mean for example if list is {1,2,3,4,5} then i have to calculate mean of Cos(1+2)+Cos(1+3)+Cos(1+4)+Cos(1+5)+Cos(2+3)+Cos(2+4)+Cos(2+5)+Cos(3+4)+Cos(3+5)+Cos(4+5)

My simple code is

     SetDirectory["/Users/macos/"]

     phia = Flatten[Import["p180.dat", "Table"]];(* for negative charges*)
     Length[phia]
     res = 0;
     l = 0;
     Do[
        Do[
          res = res + Cos[phia[[m]] + phia[[n]]];
          l++;
          , {m, 2, Length[phia]}];
          , {n, Length[phia]-1}]; // AbsoluteTiming
          mean=res/l;

it takes 5943 seconds. it is very slow. how to speed this process?[![enter image description here][2]][2] thanking in advance i also used Table command that also does not helped much. screenshot of code and time taken

$\endgroup$
  • 1
    $\begingroup$ Iterator specs to Table and Do can be done all in a single call. Your performance issue probably comes at least in part from that l++ though. You can do this much more functionally, by first building your list of phia[[n]]+phia[[m]] using Table, calling Cos on that built list (Cos is Listable), then calling Total on that. $\endgroup$ – b3m2a1 Oct 16 '17 at 15:13
16
$\begingroup$

I do not agree with @b3m2a1: functional programming does not necessarily shine here, because things involving e.g., Subsets[phia,2] require $O(n^2)$ flops and $O(n^2)$ memory, while your implementation requires $O(n^2)$ flops but only $O(n)$ memory. So, in principle, your procedural code is fine (apart from the negligible fact that it does not seem to do what you expect); the issue is that the Mathematica kernel is not good at executing procedural code without preparation. That's what Compile is good for!

So, let's compile your code. I removed the ++l, because it is not necessary.

cf = Compile[{{x, _Real, 1}}, Block[{res = 0.},
    Do[
     Do[
      If[m != n, res += Cos[x[[m]] + x[[n]]]],
      {m, 1, Length[x]}],
     {n, 1, Length[x]}];
    res/(Length[x] (Length[x] - 1))
    ],
   CompilationTarget -> "C"];

Now let's generate some data:

phia = RandomReal[{-1, 1}, 44359];

And let it run:

mean1 = cf[phia]; // AbsoluteTiming // First

(* 19.3747 *)

Still, cf requires $O(n^2)$ evaluations of trigonometric functions. We can exploit the identity $\cos(a+b) = \cos(a) \, \cos(b) - \sin(a) \, \sin(b)$ in order to cut that down to $O(n)$ trigonometric function evaluations as follows:

cf2 = Compile[{{cosx, _Real, 1}, {sinx, _Real, 1}},
   Block[{res = 0., diagres = 0.},
    Do[
     Do[
      res += (cosx[[m]] cosx[[n]] - sinx[[m]] sinx[[n]]),
      {m, 1, Length[cosx]}],
     {n, 1, Length[cosx]}
     ];
    Do[diagres += (cosx[[m]]^2 - sinx[[m]]^2), {m, 1, Length[cosx]}];
    (res - diagres)/(Length[cosx] (Length[cosx] - 1))
    ],
   CompilationTarget -> "C"];

mean2 = cf2[Cos[phia], Sin[phia]]; // AbsoluteTiming// First
mean1 == mean2

(* 2.302 *)
(* True *)

Finally, a version that also employs parallelization; it includes some rounding error, but I think it's still okay.

cf3 = Compile[{{cosx, _Real}, {sinx, _Real}, {cosy, _Real, 1}, {siny, _Real, 1}},
   cosx Total[cosy] - sinx Total[siny],
   CompilationTarget -> "C",
   RuntimeAttributes -> {Listable},
   Parallelization -> True
   ];
f3 = x \[Function] With[{cosx = Cos[x], sinx = Sin[x]},
    (Total[cf3[cosx, sinx, cosx, sinx]] - Total[cosx^2] + Total[sinx^2])/(Length[x] (Length[x] - 1))
    ];

mean3 = f3[phia]// AbsoluteTiming// First
Abs[mean2 - mean3]

(* 0.804182 *)
(* 4.05276*10^-12 *)

Running b3m2a1's implementation on this data set takes considerably longer, consumes over 50 GB(!) RAM and finally crashes the kernel. Thus, I give an example with a considerably smaller data set:

phia = RandomReal[{-1, 1}, 4435];
mean1 = cf[phia]; // AbsoluteTiming // First
mean2 = cf2[Cos[phia], Sin[phia]]; // AbsoluteTiming // First
mean3 = f3[phia]; // AbsoluteTiming // First
mean4 = Mean@Cos@Map[Total, Subsets[phia, {2}]]; // AbsoluteTiming // First

(* 0.188377 *)
(* 0.022097 *)
(* 0.008892 *)
(* 2.141 *)

Testing the accuracy:

mean2 - mean1
mean3 - mean1
mean4 - mean1

(* -2.66454*10^-15 *)
(* -7.14984*10^-14 *)
(* -5.88418*10^-14 *)

Fun fact

Due to the simplification by trigonometrics, we can obtain the same result with $O(n)$ complexity and in a completely vectorized way:

f5 = x \[Function] With[{cosx = Cos[x], sinx = Sin[x]},
    (Total[cosx]^2 - Total[sinx]^2 - Total[cosx^2] + 
       Total[sinx^2])/(Length[cosx] (Length[cosx] - 1))
    ];
mean5 = f5[phia]; // AbsoluteTiming // First
mean1 - mean5

(* 0.000195 *)
(* -1.00253*10^-13 *)

In the large example above, f5[phia] takes 0.00089 seconds; this is several million times faster than the OP's implementation.

$\endgroup$
  • $\begingroup$ You should compare the timing of your compiled version with b3m2a1's functional method (on the same machine). $\endgroup$ – bill s Oct 16 '17 at 16:36
  • $\begingroup$ Good point. Done. $\endgroup$ – Henrik Schumacher Oct 16 '17 at 17:13
  • $\begingroup$ wow. Thanks a lot Henry Schumacher, yes this is what i am doing . but when i run your code it gives me error CreateLibrary::cmperr: Compile error: LINK : fatal error LNK1104: cannot open file 'LIBCMT.lib' >> Compile::nogen: A library could not be generated from the compiled function. >> $\endgroup$ – IrfanS Oct 17 '17 at 5:01
  • $\begingroup$ Hm. I googled the error message LNK1104 and it seems to be related to Visual Studio. This sort of tells me you are running Mathematica on Windows, right? Do you have a C compiler installed on your system? Could you please run the following code on Mathematica and send me the outcome? Needs["CCompilerDriver`"]; CCompilers[] $\endgroup$ – Henrik Schumacher Oct 17 '17 at 5:45
  • 2
    $\begingroup$ You need to enclose anything with backticks in double backticks: Needs["CCompilerDriver`"]. $\endgroup$ – J. M. will be back soon Oct 25 '17 at 9:11
5
$\begingroup$

So I left a comment about this, but here's where functional programming shines. Consider this:

phia = RandomReal[{-1, 1}, 4985];
Mean@Cos@Map[Total, Subsets[phia, {2}]] // AbsoluteTiming

{3.20715, 0.70717}

If I understand correctly, that's all you want out of your code. And it's short and fast.

$\endgroup$
  • $\begingroup$ thank you . your code is also correct and no doubt i tested it is faster then my code. but need more faster you can check if you increase 4985 to 49850, because data list i use has length around ~50000 $\endgroup$ – IrfanS Oct 17 '17 at 4:55

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.