657
edits
Cookbook: Difference between revisions
Jump to navigation
Jump to search
→Distance between points
Carandraug (talk | contribs) m (→Parametrized Functions: no need to break lines when using the Math extension) |
|||
| 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 == | ||