Discussion :

[franck06400]

bonjour,

je suis actuellement stagiaire dans une société spatiale et je developpe un petit outil en matlab pour faire des calculs d'orbites, mon soucis c'est qu'une simultation prend presque 48 h !!!

j'ai en fait une fonction tres gourmande en temps de calcul qui est lancée plusieurs milliards de fois donc gratter quelques milliseconde serait un plus !

la fonction calcule un angle en radian entre 2 vecteurs (en 3 dimension.
en entrée : v et u deux vecteurs 3-sur-n
en sortie : un angle en radian entre u et v (n-sur-1)

et le code que je pense deja optimisé au max, mais si quelqu'un connait une astuce ou deux ... pour gratter du temps CPU...

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 
function angle = searchAngle(u, v)
 
norm = @(v) sqrt(sum(v.^2, 2));
dot = @(u, v) sum(u .* v, 2);
cross = @(a, b) [ a(:,2) .* b(:,3) - a(:,3) .* b(:,2), ...
      a(:,3) .* b(:,1) - a(:,1) .* b(:,3), ...
      a(:,1) .* b(:,2) - a(:,2) .* b(:,1) ];
 
normVect = norm(u) .* norm(v);
 
dotVect = dot(u, v); 
threshold = normVect * 0.9999;
 
idx1 = dotVect > threshold;
axis = cross(v(idx1,:), u(idx1,:));
angle(idx1) = asin(norm(axis) ./ normVect(idx1));
 
idx2 = dotVect < -threshold;
axis = cross(v(idx2,:), u(idx2,:));
angle(idx2) = pi - asin(norm(axis) ./ normVect(idx2));
 
idx = ~(idx1 | idx2);
angle(idx) = acos(dotVect(idx) ./ normVect(idx));
 
end

[kajtp]

Un truc qui doit sans doute ralentir ton code c'est que tu ne fais pas de "préallocation" pour la variable angle, du coup le programme doit resizer angle à la volée pour

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 
angle(idx1) = asin(norm(axis) ./ normVect(idx1)); %redimentionnement a la volee de angle, c est lent
 
angle(idx2) = pi - asin(norm(axis) ./ normVect(idx2));  %pareil

Du coup tu pourrais définir angle en début de fonction comme

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
 
angle = zeros(1,n)
 
%ensuite pareil
angle(idx1) = asin(norm(axis) ./ normVect(idx1)); % pas de resize angle est deja pret
 
angle(idx2) = pi - asin(norm(axis) ./ normVect(idx2));  %pareil

De plus dans les versiond recented de matlab les fonction cross et dot sont définies, tu devrais pouvoir les utiliser telle quelle sans les redefinir.
Sinon petit truc :

dot(a,b) = a*b'
norm(a) = sqrt(a*a')

[Kangourou]

salut,

une possibilité serait de se passer des pointeurs de fonctions. En calculant directement les vecteurs tu devrais gagner (un tout petit peu) en vitesse :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

normU = sum(u.^2, 2); normV = sum(v.^2, 2);

Pour le produit vectoriel, tu peux essayer de regarder cette fonction (un produit vectoriel qui se veut rapide, tu dois pouvoir n'utiliser que les deux dernière lignes) :
http://www.mathworks.com/matlabcentr...ectorCross3d.m

Le gros des calculs se fait de toutes facons dans asin et acos, et là c'est difficile d'accélerer plus...

Sinon il semble que tes vecteurs soient dans des tableaux N-par-3, et par 3 par N, non ?

[Jerome Briot]

Quel est l'ordre de grandeur de n ?

dot(a,b) = a*b'
norm(a) = sqrt(a*a')

Ce qui s'écrit avec MATLAB sans l'opérateur de transposition :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
a*b(:);
sqrt(a*a(:))

[franck06400]

Quel est l'ordre de grandeur de n ?

la valeur de n est 222651 !!

je regarde ton code..

[Jerome Briot]

Tu gagneras déjà du temps en utilisant des vecteurs 3xN et non pas Nx3 :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
clc
clear all
 
n = 222651;
 
u = rand(n, 3);
v = rand(n, 3);
 
tic
angle = searchAngle(u, v);
toc
 
clear functions
clear angle
 
tic
u = u.';
v = v.';
angle = searchAngle(u, v);
toc

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
Elapsed time is 0.026036 seconds.
Elapsed time is 0.018571 seconds.

L'opérateur de transposition est juste la pour l'exemple.
L'idéal serait de pouvoir modifier les dimensions des deux tableaux dès leur création.