Prev: a question about angle in solve
Next: Bloated files when saving large variables to .mat files with -append
From: Torsten Hennig on 6 May 2008 22:48
>Dear all,
>
>I am not an expert in optimization and hence I'd really
>appreciate your help with following problem.
>
>Consider n random sampling locations in 2-D. E.g.
>
>n = 100;
>x = rand(n,2);
>
>The pairwise euclidean distance between all sampling >points
>is (in case you have the Stats toolbox)
>
>d = pdist(x);
>
>You may want to plot a histogram of d
>
>hist(d,20);
>
>and realize, that the distances have a non-uniform
>distribution. But a uniform distribution of the sampling
>points is what I want. Hence, I'd like to shift the >sampling
>locations so that their pairwise distance approximates a
>continuous uniform distribution.
>
>I have no idea how to formulate such an optimization
>problem. If you have an idea, hint, ... whatever, >that'll be
>great. If the question is too general to be answered >here,
>please let me know. This is no homework :-)
>
>Best regards,
>Wolfgang


Hi Wolfgang,

maybe you should tackle the problem 'the other way round'.
Given that the points are to be distributed within the
rectangle x_min <= x <= x_max, y_min <= y <= y_max,
you could first generate distances d_ij (1<=i<=n,
i+1<=j<=n) uniformly distributed over
[0; sqrt((x_max-x_min)^2 + (y_max-y_min)^2)].
After this, you could calculate a series of
points (x_i,y_i) (1<=i<=n) by solving the following
optimization problem:
minimize : eps
-eps <= d_ij^2 - (x_i-x_j)^2 - (y_i-y_j)^2 <= eps
(1<=i<=n, i+1<=j<=n)
x_min <= x_i <= x_max (1<=i<=n)
y_min <= y_i <= y_max (1<=i<=n)

If a solution with a 'small' value of eps exists, you
can be sure that the pointwise distances are
distributed quite uniformly.

It's an interesting problem ; it's worth a try to
post it to sci.math as well.

Best wishes
Torsten.
From: Wolfgang Schwanghart on 7 May 2008 09:24
Thanks a lot to you all for your helpful comments. I wrote a
function for the problem as suggested by Torsten. So far
it's rarely an optimization of a 2-D geostatistical sampling
scheme, because nearly always I get a circle type of scheme.
And also it might be most economical to cover all ranges of
distance within a sampling domain, I doubt its acceptance
among soil scientists. :-)

Best regards,
Wolfgang




function [x,x0,dnow,d] = samplescheme(n,dom)

% optimization of 2-D geostatistical sampling scheme
%
% [x,x0,dnow,d] = samplescheme(n,dom)
% ___________________________________________________
%
% Generate a point sampling scheme where the
% distribution of the interpoint distances follows
% a theoretical distribution.
%
% Input:
% n: number of desired sampling locations
% (scalar)
% dom: boundaries of rectangular sampling
% domain ([[minx maxx miny maxy])
%
% Output:
% x: sampling locations (nx2 matrix)
% x0: initial random configuration of sampling
% locations (nx2 matrix)
% dnow: vector with pairwise distances
% d: the desired pairwise distances
%
% Requirements:
% Optimization Toolbox, Statistical Toolbox
%
% Example:
%
% [x,x0,dnow,d] = samplescheme(50,[0 1 0 1])
% scatter(x0(:,1),x0(:,2),'.')
% hold on
% scatter(x(:,1),x(:,2),'rs')
% axis image
% legend('x0','x')
% figure
% subplot(1,2,1);
% hist(d,20);
% title('theor. distribution of sample spacings');
% subplot(1,2,2);
% hist(dnow,20);
% title('emp. distribution of sample spacings');
%
% ___________________________________________________
%



% nr of pairwise distances
n_d = (n^2-n)/2;

% upper and lower limits of distances
dlu = [0.001; ...
sqrt((dom(2)-dom(1))^2 + (dom(4)-dom(3))^2)];

% create desired distance distribution here
% d = rand(n_d,1);
% d = raylrnd(ones(n_d,1));
% d = log(linspace(1,10,n_d)');
d = hypot(rand(n_d,1),rand(n_d,1));

% linear transformation of d to fit upper and lower
% domain boundaries
d = (d-min(d))*(dlu(2)-dlu(1))/max(d-min(d))+dlu(1);

% indices for pairwise distances
IX1 = 1:n;
IX2 = (1:n)';

[IX1,IX2] = meshgrid(IX1,IX2);
[ignore1,ignore2,IX1] = find(spdiags(IX1, -n:-1));
[ignore1,ignore2,IX2] = find(spdiags(IX2, -n:-1));

clear ignore1 ignore2

ij = [IX1(:) IX2(:)];

clear IX1 IX2

% initial values
x0 = rand(n,2);

x0(:,1) = x0(:,1) * (dom(2)-dom(1)) + dom(1);
x0(:,2) = x0(:,2) * (dom(4)-dom(3)) + dom(3);

lb = [repmat(dom(1),n,1) repmat(dom(3),n,1)];
ub = [repmat(dom(2),n,1) repmat(dom(4),n,1)];

% function handle
fun = @(x) d.^2 - ...
(x(ij(:,1),1)-x(ij(:,2),1)).^2 - ...
(x(ij(:,1),2)-x(ij(:,2),2)).^2;

% use lsqnonlin to solve the problem
x = lsqnonlin(fun,x0,lb,ub);
dnow = pdist(x)';






First  |  Prev  | 
Pages: 1 2 3
Prev: a question about angle in solve
Next: Bloated files when saving large variables to .mat files with -append