function [datamatrix_cb,datamatrix_eu,corr_cb_eu,gamma_cb_eu,datavectors,coord] ...
    = data_generate_cityblock_and_euclidean(n,error,elongation_factor)

% DATA_GENERATE_CITYBLOCK_AND_EUCLIDEAN produces data having an underlying
% city-block and Euclidean structure.
%
% syntax: [datamatrix_cb,datamatrix_eu,corr_cb_eu,gamma_cb_eu,datavectors,coord] ...
%    = data_generate_cityblock_and_euclidean(n,error,elongation_factor)
%
% The input parameters are the number of objects, $n$; ERROR is in the form of
% a multiplier used on a randomly (entry) permutated proximity matrix with
% perfect city-block or Euclidean structure that is then added to these
% latter matrices; the ELONGATION_FACTOR refers to a multiplier used on the
% randomly generated first dimension coordinates (in case a dominant
% dimension is desired).
% The outputs are the cityblock and Euclidean $n \times n$ proximity
% matrices, DATAMATRIX_CB and DATAMATRIX_EU, respectively; CORR_CB_EU and
% GAMMA_CB_EU are the Pearson and Goodman-Kruskal coefficients of correlation 
% between the entries in DATAMATRIX_CB and DATAMATRIX_EU; DATAVECTORS is
% $n*(n-1)/2 \times 2$ and contains the entries in DATAMATRIX_CB and
% DATAMATRIX_EU in vector form; COORD contains the randomly generated
% coordinates.


rand('state', sum(100*clock));

coord = rand(n,2);
coord(:,1) = coord(:,1)*elongation_factor;
datavectors = zeros(n*(n-1)/2,2);
datamatrix_cb = zeros(n,n);
datamatrix_eu = zeros(n,n);

for i = 1:(n-1)
    for j = (i+1):n
        
        datamatrix_cb(i,j) = abs(coord(i,1) - coord(j,1)) + ...
            abs(coord(i,2) - coord(j,2));
        
        datamatrix_eu(i,j) = sqrt(((coord(i,1) - coord(j,1))^2) + ...
            ((coord(i,2) - coord(j,2))^2));
        
        datamatrix_cb(j,i) = datamatrix_cb(i,j);
        datamatrix_eu(j,i) = datamatrix_eu(i,j);
        
    end
end

[randprox_cb] = proxrand(datamatrix_cb);
[randprox_eu] = proxrand(datamatrix_eu);

datamatrix_cb = datamatrix_cb + (error*randprox_cb);
datamatrix_eu = datamatrix_eu + (error*randprox_eu);

index = 0;


for i = 1:(n-1)
    for j = (i+1):n
        
        index = index + 1;
        
        datavectors(index,1) = datamatrix_cb(i,j);
        datavectors(index,2) = datamatrix_eu(i,j);
        
      
        
    end
end



corr_cb_eu = corrcoef(datavectors);

length_vector = n*(n-1)/2;
splus = 0.0;
sminus = 0.0;

for i = 1:(length_vector - 1)
    for j = (i+1):length_vector
        
        signfunction = (datavectors(i,1) - datavectors(j,1))* ...
            (datavectors(i,2) - datavectors(j,2));
        
        if(signfunction > 0.0)
            
            splus = splus + 1;
            
        end
        
        if(signfunction < 0.0)
            
            sminus = sminus + 1;
            
        end
    end
end

gamma_cb_eu = (splus - sminus)/(splus+sminus);

figure(1)

plot(datavectors(:,1),datavectors(:,2),'b.')


