# Efficiency in vector array manipulation

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.

• Can you describe what you are trying to accomplish here? Jul 18 '17 at 15:27
• 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. Jul 18 '17 at 15:47
• 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]. Jul 18 '17 at 18:18

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 //DeveloperPackedArrayQ


{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 = DeveloperToPackedArray @ 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}

• Wow, thank you so much @CarlWoll. A single plot when from 15 + min to 2.7 seconds. Huge improvement! Jul 18 '17 at 18:03
• 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? Jul 18 '17 at 18:06
• Probably your data has a mix of integers and reals. You could always try something like ndat = DeveloperToPackedArray @ N @ dat and then see if DeveloperPackedArrayQ @ ndat is true. Jul 18 '17 at 18:11
• @CarlWollThat solved the problem, Thank you. Jul 18 '17 at 18:22
• 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) Jul 18 '17 at 18:26