--- solve_problems.m
 
function report = solve_problems()
% data_files={'Alberta';'Galvao100';'Galvao100';'Galvao150';'Galvao150';'Galvao150'};
% p=[10 10 15 5 15 20];
% data_files={'TestData'};
% p=[2];
data_files={'Alberta'};
p=[10];
for i=1:length(data_files)
    data_file=char(data_files(i));
    dist=distances([data_file '_distances.txt']);
    demand=load([data_file '_demands.txt']);
    [bestLB,iterations,debug]=solve_p_median(dist,demand,p(i));
    report(i).bestLB=bestLB;
    report(i).iterations=iterations;
    report(i).debug=debug;
    dlmwrite([data_file '_p' int2str(p(i)) '_iterations.txt'],iterations);
    dlmwrite([data_file '_p' int2str(p(i)) '_debug.txt'],debug);
end
 
--- distances.m
 
function distances=distances(data_file)
distance_data=load(data_file);
distances=zeros(distance_data(end,1)+1);
for row=1:length(distance_data)
    from=distance_data(row,1);
    to=distance_data(row,2);
    distance=distance_data(row,3);
    distances(to,from)=distance;
    distances(from,to)=distance;
end
 
--- solve_p_median.m
 
function [bestLB,iterations,debug]=solve_p_median(dist,demand,p)
bestLB=0;
bestUB=inf;
currentLB=0;
currentUB=inf;
iterations(1,:)=[0 currentLB bestLB currentUB bestUB];
pi=2;
n_c=length(demand); % number of customers
n_s=length(dist(:,1));  % number of sites
u=ones(n_c,1); % LR (Lagrangean Relaxation) multipliers
zero=zeros(n_c,1);
x=zeros(n_s,n_c); % site/customer assignments
i=1;
piUpdateTime=1;
improvementOccurred=0;
debug=[0 0 2 zeros(1,p) u']; % parameters stored for debugging (iteration, step_size, pi, open facilities, u)
while(~stoppingCondition(pi,bestLB,bestUB,i))
    for s=1:n_s
        cost=dist(:,s).*demand;
        newCost=cost-u;
        z_LR(s)=sum(min(zero,newCost));
        x(s,find(newCost<0))=1; % x_{s,c} = 1 if cost negative
    end
    [z_sorted,order]=sort(z_LR);
    currentLB=sum(z_sorted(1:p))+sum(u); % z_{LR}(u) value is current LB
    facilities=order(1:p); % open p facilities where z is smallest
    x(setdiff(1:n_s,facilities),:)=0; % x_{s,c} \leq y_s \forall s,c constraint
    if(currentLB>bestLB)
        bestLB=currentLB;
    end
    currentUB=findUB(facilities,dist,demand);
    if(currentUB<bestUB)
        bestUB=currentUB;
    end
    iterations(i+1,:)=[i currentLB bestLB currentUB bestUB];
    normOfRelaxedCsts=sum((1-sum(x)).^2);
    if(normOfRelaxedCsts == 0) % hit the lower bound 
        break
    end
    step=pi*(bestUB-currentLB)/normOfRelaxedCsts; % s^t = {\pi (UB* - z_{LR}(u^t)) \over \sum_c (1-\sum_s x_{sc})^2}
    u=u+step*(1-sum(x))'; % u_c^{t+1} = u_c^t + s^t (1-\sum_s x_{sc}^t)
    if(~improvementsOccur(iterations,piUpdateTime))
        pi=pi/2;
        piUpdateTime=i;
    end
    debug(i+1,:)=[i step pi facilities u'];
    i=i+1
end
 
function result=stoppingCondition(pi,bestLB,bestUB,iterationNo)
result = (pi < 0.005) | (bestUB == bestLB);
% result = iterationNo > 15000 | (bestUB == bestLB);
 
function result=improvementsOccur(iterations,piUpdateTime)
n=30;
currentTime=length(iterations(:,1));
timeSinceLastPiUpdate=currentTime-piUpdateTime;
if(currentTime <= n | timeSinceLastPiUpdate <= n)
    result = 1;
else
    lastImpForLB=whenDidLastImprovementOccur(iterations(end-n:end,3));
    lastImpForUB=whenDidLastImprovementOccur(iterations(end-n:end,5));
    if(lastImpForLB <= n | lastImpForUB <= n)
        result = 1;
    else
        result = 0;
    end
end
 
function lastImp=whenDidLastImprovementOccur(iterations)
lastValue=iterations(end);
ixLastImp=length(find(iterations~=lastValue));
lastImp=length(iterations)-ixLastImp;
 
function currentUB=findUB(facilities,dist,demand)
feasibleAssignments=assignCustomers(facilities,dist);
currentUB=sum(feasibleAssignments.*dist)*demand;
 
function customerAssignments=assignCustomers(facilities,dist)
n_s=length(dist(:,1)); % number of sites
n_c=length(dist(1,:)); % number of customers
customerAssignments=zeros(n_s,n_c);
distOpenFacilities=dist(facilities,:);
[minDist,order]=min(distOpenFacilities);
for c=1:n_c
    customerAssignments(facilities(order(c)),c)=1;
end