2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Henrik Schumacher Nov 24 '17 at 17:39
3
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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!

$\endgroup$

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.