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

% CITYBLOCK_MULTIFACILITY_HAP carries out the Cityblock multifacility
% location problem using a hyperboloid approximation procedure.
%
% syntax: [new_facility_location,distance_new_to_old,...
%    distance_among_new,objective_function] = ...
%    cityblock_multifacility_hap(weight_new_to_old,...
%    weight_among_new,old_facility_location,epsil)
%
% 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; and EPSIL (e.g., 1.0e-012) produces the hyperboloid
% approximation allowing the objective function to be differentiable
% everywhere.
% 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);

new_facility_location = zeros(numdim,numnew);
weightsum = sum(weight_new_to_old');

for d = 1:numdim
    for j = 1:numnew
        
        if (weightsum(j) > 0)
            
            for i = 1:numold
                
                new_facility_location(d,j) = new_facility_location(d,j) + ...
                (weight_new_to_old(j,i)*old_facility_location(d,i))/weightsum(j);
                
            end
            
        else
            
            for i = 1:numold
                
                new_facility_location(d,j) = new_facility_location(d,j) + ...
                old_facility_location(d,i)/numold;
                
            end
        end
    end
end

new_facility_location;


distance_among_new = zeros(numnew,numnew);
distance_new_to_old = zeros(numnew,numold);
dim_distance_among_new = zeros(numdim,numnew,numnew);
dim_distance_new_to_old = zeros(numdim,numnew,numold);

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

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

objective_function = 0;

for j = 1:(numnew-1)
    for k = (j+1):numnew
        
        for d = 1:numdim
        
        objective_function = objective_function + weight_among_new(j,k)* ...
            sqrt((dim_distance_among_new(d,j,k)^2 + epsil));
        
        end
        
    end
end

for j = 1:numnew
    for k = 1:numold
        
        for d = 1:numdim
        
        objective_function = objective_function + weight_new_to_old(j,k)* ...
            sqrt((dim_distance_new_to_old(d,j,k)^2 + epsil));
        
        end
        
    end
end

obdiff = 1.0;

iterations = 0;

while (obdiff > 1.0e-009)
    
    iterations = iterations + 1;
    
    obprev = objective_function;
    
    new_facility_location_old = new_facility_location;
    
    new_facility_location = zeros(numdim,numnew);
    
    for d = 1:numdim
        for j = 1:numnew
            
            for i = 1:numold
                
                new_facility_location(d,j) = new_facility_location(d,j) + ...
                    (weight_new_to_old(j,i)*old_facility_location(d,i))/ ...
                    (sqrt((dim_distance_new_to_old(d,j,i)^2 + epsil)));
                
            end
            
            for k = 1:numnew
                
                if (j ~= k)
                    
                    new_facility_location(d,j) = new_facility_location(d,j) + ...
                     (weight_among_new(j,k)*new_facility_location_old(d,k))/ ...
                     (sqrt((dim_distance_among_new(d,j,k)^2 + epsil)));
                 
                end
            end
        end
    end
    
    dim_weightsum = zeros(numnew,numdim);
    
    for j = 1:numnew
        for i = 1:numold
            
            for d = 1:numdim
            
            dim_weightsum(j,d) = dim_weightsum(j,d) + weight_new_to_old(j,i)/ ...
                (sqrt((dim_distance_new_to_old(d,j,i)^2 + epsil)));
            
            end
        end
        
        for k = 1:numnew
            
            if (j ~= k)
                
                for d = 1:numdim
                
                dim_weightsum(j,d) = dim_weightsum(j,d) + weight_among_new(j,k)/ ...
                    (sqrt((dim_distance_among_new(d,j,k)^2 + epsil)));
                
                end
            end
        end
    end
    
    for d = 1:numdim
        for j = 1:numnew
            
            new_facility_location(d,j) = new_facility_location(d,j)/dim_weightsum(j,d);
            
        end
    end
    
    
   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
            
            dim_distance_among_new(d,j,k) =  abs(new_facility_location(d,j) - new_facility_location(d,k));
            
            distance_among_new(j,k) = distance_among_new(j,k) + ...
            abs(new_facility_location(d,j) - new_facility_location(d,k));
            
        end
        
      
        
    end
end

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

objective_function = 0;

for j = 1:(numnew-1)
    for k = (j+1):numnew
        
        for d = 1:numdim
        
        objective_function = objective_function + weight_among_new(j,k)* ...
            sqrt((dim_distance_among_new(d,j,k)^2 + epsil));
        end
        
    end
end

for j = 1:numnew
    for k = 1:numold
        
        for d = 1:numdim
        
        objective_function = objective_function + weight_new_to_old(j,k)* ...
            sqrt((dim_distance_new_to_old(d,j,k)^2 + epsil));
        
        end
        
    end
end
    
    obdiff = abs(obprev - objective_function)/obprev;
    
end

iterations;

%toc;
        
    
        
