# Changes

,  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 [itex] (x_i,y_i,\ldots)[/itex]. 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