2
$\begingroup$

I have a plotting subroutine that contains the following two steps. Im pretty sure this is not he most efficient implementation. The problematic steps are the following:

1)

dat = RandomReal[{0, 2}, {400000, 3}];
k=Dimension[dat][[1]]; 
z0=2.5;
Parallelize[list = Reap[Do[
    Sow[{dat[[n, 1]], dat[[n, 2]], dat[[n, 3]]/z0}];
    Sow[{dat[[n, 1]], dat[[n, 2]] - 1.154701, dat[[n, 3]]/z0}];
    Sow[{dat[[n, 1]] - 1., dat[[n, 2]] - 0.577350, 
      dat[[n, 3]]/z0}];
    Sow[{dat[[n, 1]] - 1., dat[[n, 2]] - 1.154701 - 0.577350, 
      dat[[n, 3]]/z0}];
    , {n, 1, k}]
   ][[2]][[1]];]; 

2)

regularPolygon[nbrSides_Integer?(# > 2 &), scale_:1]:= 
Polygon[scale{Sin[#], -Cos[#]} & /@ (2 Pi*Range[1/nbrSides, 1, 1/nbrSides])];

listhex = Reap[Do[If[
   Element[{list[[m, 1]], list[[m, 2]]},regularPolygon[6, 0.645497]] ==True,Sow[list[[m]]], 
   Continue], {m, 1, 4*k}]][[2]][[1]];

With the current dimensions of the array these step take a very long time.

$\endgroup$
  • $\begingroup$ Can you describe what you are trying to accomplish here? $\endgroup$ – MarcoB Jul 18 '17 at 15:27
  • $\begingroup$ Sure, so I have a set of vectors, the original data,I then make a larger list by doing 3 translation to the original set.I then use the regularPolygon to select the points inside a hexagon and then do a ListDensityPlot of that data at the end. $\endgroup$ – Giovanni Baez Jul 18 '17 at 15:47
  • 2
    $\begingroup$ Your regularPolygon is equivalent to the built-in Polygon@CirclePoints[{scale, 3 Pi/2 + 2 Pi/nbrSides}, nbrSides], or effectively equivalent to the simpler Polygon@CirclePoints[{scale, 3 Pi/2}, nbrSides]. Or as a region, it's equivalent to RegularPolygon[{scale, 3 Pi/2}, nbrSides]. $\endgroup$ – Michael E2 Jul 18 '17 at 18:18
5
$\begingroup$

Using Parallelize to do basic matrix operations is a terrible idea, because Mathematica already uses very effective parallelization under the hood for such operations. Here is an alternate approach for creating your list object. First, your data:

z0 = 2.5;
dat = RandomReal[{0, 2}, {400000, 3}];

Next, I would prepend each column with a 1. so that I can use Dot to construct your new lists:

datOne = Join[
    ConstantArray[1., {400000, 1}],
    dat,
    2
]; //AbsoluteTiming

datOne //Developer`PackedArrayQ

{0.00737, Null}

True

Note that the new matrix datOne is still packed. Now, here is a transformation matrix that will take an element {a, b, c} of dat and convert it into {a - 1., b - .57735, c/z0}:

tm3 = Developer`ToPackedArray @ N @ {{-1, -.57735, 0}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1/z0}};

For example:

{{1,a,b,c}} . tm3

{{-1. + 1. a, -0.57735 + 1. b, 0. + 0.4 c}}

So, to generate the 3rd Sow of your step 1, you would just do:

s3 = datOne . tm3; //AbsoluteTiming

{0.004051, Null}

The other sows can be handled in the same way, and in my answer I will ignore them. You could just repeat this procedure 4 times with 4 different transformation matrices.

Now, for step 2. Instead of using Element[point, region] to test for region membership, it will be an order of magnitude faster to construct a RegionFunction, and then apply that RegionFunction to all of your points: For example:

regularPolygon[nbrSides_Integer?(# > 2 &), scale_:1] := Polygon[
    scale{Sin[#], -Cos[#]} & /@ (2 Pi*Range[1/nbrSides, 1, 1/nbrSides])
];

rf = RegionMember[regularPolygon[6, .645497]]; //AbsoluteTiming

{0.06491, Null}

Applying the RegionFunction to your data (notice that rf is not being mapped over the data points):

bool = rf[s3[[All, ;;2]]]; //AbsoluteTiming
Tally[bool]

{0.025368, Null}

{{False, 292548}, {True, 107452}}

To obtain the list of in-hexagon points you can use Pick:

listhex = Pick[s3, bool]; //AbsoluteTiming

{0.189166, Null}

| improve this answer | |
$\endgroup$
  • $\begingroup$ Wow, thank you so much @CarlWoll. A single plot when from 15 + min to 2.7 seconds. Huge improvement! $\endgroup$ – Giovanni Baez Jul 18 '17 at 18:03
  • $\begingroup$ I did notice something when using my real data, even though is in the same format as dat, {{x1,y1,z1},{x2,y2,z2},...}, it returns false when it checks for packedarray. I import my data from a text file and then make the array in the format I described. Do you know what is the issue? $\endgroup$ – Giovanni Baez Jul 18 '17 at 18:06
  • 1
    $\begingroup$ Probably your data has a mix of integers and reals. You could always try something like ndat = Developer`ToPackedArray @ N @ dat and then see if Developer`PackedArrayQ @ ndat is true. $\endgroup$ – Carl Woll Jul 18 '17 at 18:11
  • $\begingroup$ @CarlWollThat solved the problem, Thank you. $\endgroup$ – Giovanni Baez Jul 18 '17 at 18:22
  • 4
    $\begingroup$ FWIW, datOne = PadLeft[dat, {Automatic, 4}, 1.] is about 5-10% faster than Join + ConstantArray -- practically negligible in this case though. But it's also shorter to type and you can avoid having to input Length@dat. (+1) $\endgroup$ – Michael E2 Jul 18 '17 at 18:26

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.