function [new_facility_location,distance_new_to_old,...
    distance_among_new,objective_function] = ...
    squared_euclidean_multifacility(weight_new_to_old,...
    weight_among_new,old_facility_location)

% SQUARED_EUCLIDEAN_MULTIFACILITY carries out the squared Euclidean multifacility
% location problem using closed form expressions for the results.
%
% syntax: [new_facility_location,distance_new_to_old,...
%    distance_among_new,objective_function] = ...
%    squared_euclidean_multifacility(weight_new_to_old,...
%    weight_among_new,old_facility_location)
%
% There are n new facilities to place among m existing facilities in d
% dimensions.  WEIGHT_NEW_TO_OLD (n x m) and WEIGHT_AMONG_NEW (n x n) are
% both weight matrices (with similarity interpretations);
% OLD_FACILITY_LOCATION (d x m) contains the coordinates for the existing
% facilities. 
% As output, there are the NEW_FACILITY_LOCATIONS (d x n);
% DISTANCE_NEW_TO_OLD (n x m); DISTANCE_AMONG_NEW (n x n); and the
% OBJECTIVE_FUNCTION (i.e., the weighted sum of the distances among the new
% locations plus of the new locations to those existing).


%tic;

numnew = size(weight_new_to_old,1);
numold = size(weight_new_to_old,2);
numdim = size(old_facility_location,1);

weighted_matrix = zeros(numnew,numnew);

sum_weight_new_to_old = sum((weight_new_to_old)');
sum_weight_among_new = sum(weight_among_new);

for j = 1:numnew
    for k = 1:numnew
        
        if( j ~= k )
            
            weighted_matrix(j,k) = (-1.0)*weight_among_new(j,k);
            
        else
            
            weighted_matrix(j,j) = sum_weight_new_to_old(j) + ...
                sum_weight_among_new(j);
        end
    end
end


new_facility_location = ((inv(weighted_matrix))*(weight_new_to_old)*((old_facility_location)'))';

distance_among_new = zeros(numnew,numnew);
distance_new_to_old = zeros(numnew,numold);

for j = 1:numnew
    for k = 1:numnew
        
        for d = 1:numdim
            
            distance_among_new(j,k) = distance_among_new(j,k) + ...
            (new_facility_location(d,j) - new_facility_location(d,k))^2;
            
        end
        
        
    end
end

for j = 1:numnew
    for k = 1:numold
        
        for d = 1:numdim
            
            distance_new_to_old(j,k) = distance_new_to_old(j,k) + ...
            (new_facility_location(d,j) - old_facility_location(d,k))^2;
            
        end
       
        
    end
end

objective_function = 0;

for j = 1:(numnew-1)
    for k = (j+1):numnew
        
        objective_function = objective_function + weight_among_new(j,k)* ...
            distance_among_new(j,k);
        
    end
end

for j = 1:numnew
    for k = 1:numold
        
        objective_function = objective_function + weight_new_to_old(j,k)* ...
            distance_new_to_old(j,k);
        
    end
end

%toc;
        
    
        
