Changes

Jump to navigation Jump to search
2,280 bytes added ,  15:07, 22 August 2012
Line 98: Line 98:  
=== Distance between points ===
 
=== Distance between points ===
 
==== Problem  ====
 
==== Problem  ====
 +
Given a set of points in space we want to calculate the distance between all of them. Each point is described by its components <math> (x_i,y_i,\ldots)</math>. Asusme that the points are saved in a matrix '''<tt>P</tt>''' with '''<tt>N</tt>''' rows (one for each point) and '''<tt>D</tt>''' columns, one for each component.
 
==== Solution ====
 
==== Solution ====
 +
One way of proceeding is to use the broadcast properties of operators in GNU Octave. The square distance between the points can be calculated with the code
 +
<!-- {{SyntaxHighlight| -->
 +
{{Code|Calculate square distance between points|<syntaxhighlight lang="octave" style="font-size:13px">
 +
[N, dim] = size (P);
 +
Dsq    = zeros (N);
 +
for i = 1:dim
 +
  Dsq += (P(:,i) - P(:,i)').^2;
 +
endfor
 +
</syntaxhighlight>}}
 +
This matrix is symmetric with zero diagonal.
 +
 +
Similarly the vectors pointing from one point to the another is
 +
<!-- {{SyntaxHighlight| -->
 +
{{Code|Calculate radius vector between points|<syntaxhighlight lang="octave" style="font-size:13px">
 +
R    = zeros (N,N,dim);
 +
for i = 1:dim
 +
  R(:,:,i) = P(:,i) - P(:,i)';
 +
endfor
 +
</syntaxhighlight>}}
 +
 +
The relation between <tt>Dsq</tt> and <tt>R</tt> is
 +
<!-- {{SyntaxHighlight| -->
 +
{{Code|Calculate square distance from radius vector between points|<syntaxhighlight lang="octave" style="font-size:13px">
 +
Dsq = sumsq (R,3);
 +
</syntaxhighlight>}}
 +
 +
 
==== Discussion ====
 
==== Discussion ====
 +
The calculation can be implemented using functions like <tt>cellfun</tt> and avoid the loop over components of the points. However in most cases we will have more points than components and the improvement, if any, will be minimal.
 +
 +
Another observation is that the matrix Dsq is symmetric and we could store only the lower or upper triangular part. To use this optimization in a practical way check the help of the functions <tt>vech</tt> and <tt>unvech</tt> (this one is in the Forge package ''general''). Two functions that haven't seen the light yet are <tt>sub2ind_tril</tt> and <tt>ind2sub_tril</tt> (currently private functions in the [[Mechanics_package | Forge package mechanics]]) that are useful to index the elements of a vector constructed with the function <tt>vech</tt>. Each page (the third index) of the multidimensional array <tt>R</tt> is an anti-symmetric matrix and we could also save some memory by keeping only one of the triangular submatrices.
    
== Plotting ==
 
== Plotting ==
 
== User input ==
 
== User input ==
657

edits

Navigation menu