# Vectortransformation inside Cases condition very slow

I have the following piece of code

all[phi_, theta_, alpha_] :=
EulerMatrix[{-(phi +
90) Degree, theta Degree, -(180 - alpha) Degree}, {3, 2, 1}]

condition = Cases[data, { _, _, x_, _, _, vx_, vy_, vz_, _} /;
85 < VectorAngle[all[Phi[x], Theta[x], Alpha[x]].{0, 0, -1},
{vx, vy, vz}]/Degree < 95]


with Phi[x], Theta[x] and Alpha[x] defined as Piecewise functions that give back angles in degrees along x. "data" is a list of particles, each with different positions, velocities and other properties.

What the code does is filter out all particles in data, whose velocity vector is out of a certain field of view.

Since the orientation of the field of view is not constant, but changing depending on the x coordinate inside my reference system, I need to apply a transformation matrix to the field of view-system inside my condition.

VectorAngle[all[Phi[x], Theta[x], Alpha[x]].{0, 0, -1},{vx, vy, vz}]


But the way I'm doing it seems very slow. For a 500k long list, it takes about 10 min. The main problem seems to be the calling of the rotation matrix inside VectorAngle. Is there a faster way?

Edit:

An example of the entries in data :

{0, 3.5104e-07, 45.351, -0.654687, 0.300059, 6.81731, -1.25439, -4.13684, 8.07234}

The first entry identifies the particles species, the second is a numerical weight, the last entry is actually the magnitude of the particles velocity.

Here's how I defined Phi[x] (I didn't include all conditions):

Phi[x_] :=
Piecewise[{{140, x < -14.},{125.,-14. <= x < -12.},{90.,12.<= x <-9.5},
{160., -9.5 <= x < -2.},{180.,-2. <= x < 4.5}}


Alpha[x]and Theta[x]are defined in the exact same way, just other numbers. And the intervals covered in the conditions do not overlap completely inbetween the three.

• Can you please provide a snippet of data in order to give an impression how this looks like? What are the first two and the last entry? And what are Phi[x], Theta[x], and Alpha[x]. It would be nice to have their definitions. Nov 24, 2017 at 17:39

Just some first guesses:

First, VectorAngle is notoriously slow; that's something that actually rather inacceptable. Here is a compiled version

cAngle = Compile[{{u, _Real, 1}, {v, _Real, 1}},
Block[{lu, lv},
lu = Sqrt[Abs[u.u]];
lv = Sqrt[Abs[v.v]];
If[lu > 1. 10^-12 && lv > 1. 10^-12,
ArcCos[u.v/lu/lv],
0.5 Pi
]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
]


Second, Degree is nice to have but I it may also slow things down. So I removed it.

Third, the code about the EulerMatrix needs only to be applied once symbolically. (Quite likely, it is the slowest part in the code). I put this together with a lot of other things to be done into a CompiledFunction. If Phi, Theta, and Angle should be really only Piecewise, then this should work well with Compile. (Not having your code, I cannot check this... too bad.)

Phi = x \[Function] x;
Theta = x \[Function] x;
Alpha = x \[Function] x;
With[{code =
Simplify[
EulerMatrix[{-(Phi[x] + 90) Degree,
Theta[x] Degree, -(180 - Alpha[x]) Degree}, {3, 2, 1}
].{0, 0, -1} /. Degree -> Pi/180.]},
u = Compile[{{x, _Real}},
code,
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
]
];


Some random example (which is also pretty meaningless without real example data... again too bad). I use Clip to set all entries outside of the desired interval to -42. (all other entries remain unchanged). Afterwards, Cases can detect detect the entries different from -42.

data = RandomReal[{-Pi, Pi}, {5000000, 9}];
condition = Cases[
Clip[cAngle3[u[data[[All, 3]]], data[[All, 6 ;; 8]]], {85 Pi/180.,
95 Pi/180.}, {-42., -42.}],
Except[-42.]
]; // AbsoluteTiming


{1.08703, Null}

Runs in roughly a second. Unfortunately, I cannot check if this is the correct result. Guess why. No example data.

Sadly I can't comment my post because I made it before registering. The answer of @Henrik Schumacher works very well indeed and should be approved. Here's how I applied his ideas:

condition =
Cases[data, { _, _, x_, _, _, vx_, vy_, vz_, _} /;
85 < cAngle[u[x], {vx, vy, vz}]/Degree < 95];


using his definitions for uand cAngle. This gives me exactly what I wanted, a list of the remaining elements in data` that fulfill the asked condition, done in negligible time (about 5 seconds). Thank you very much!