function [report,result] = solve_problems()
data_files={'bon287';'p654'};
p=[2 4 5 6];
n_runs=5;
k=0;
for run=1:n_runs
   for i=1:length(data_files)
      for j=1:length(p)
          k=k+1;
          data_file=char(data_files(i));
          locations=load([data_file '_locations.txt']);
          demands=load([data_file '_demands.txt']);
          [z,x,cycles,debug]=solve_location_allocation(locations,demands,p(j));
          dlmwrite(['results\' data_file '_p' int2str(p(j)) '_run' int2str(run) '_cycles.txt'],cycles);
          dlmwrite(['results\' data_file '_p' int2str(p(j)) '_run' int2str(run) '_debug.txt'],debug);
          dlmwrite(['results\' data_file '_p' int2str(p(j)) '_run' int2str(run) '_result.txt'],[z x]);
          result(run,i,j)=z;
       end
   end
end
 
function [z_best,x_best,cycles,debug]=solve_location_allocation(custLocs,demands,p,acceptanceRule)
if nargin < 4
    acceptanceRule=@probabilisticAcceptance;
end
n_c=length(demands); % number of customers
tc=0; % termination counter
ip=0.3; % updates per cycle bound for temperature
r=0.9; % cooling ratio
fp=0.02; % updates per cycle bound for termination
gamma=1.1; % cycle growth factor
p_1=0.95; % initial probability
k=1; % number of changes in the neighborhood function
thresholdParameter=0.1;
T=initializeTemperature(p_1,p,k,custLocs,demands);
ns=nchoosek(n_c,k)*(p-1)^k;
L=2*ns;
x=generateArbitrarySolution(n_c,p); % customer/facility assignments
z_best=inf;
cycle=1;
cycles(1,:)=[z_best L T tc]; % parameters stored for debugging per cycle
debug(1,:)=[z_best inf inf 0 L 0]; % parameters stored for debugging (per iteration)
while(~stoppingCondition(tc,cycles))
    j=0;
    for i=1:L
        facLocs=single_facility_optimization(x,custLocs,demands);
        x=findAssignments(facLocs,custLocs);
        z=f(x,demands,facLocs,custLocs);
        x_=pickANeighbor(x,k,p);
        facLocs_=single_facility_optimization(x_,custLocs,demands);
        x_=findAssignments(facLocs_,custLocs);
        z_=f(x_,demands,facLocs_,custLocs);
        delta=z-z_;
        if(delta<0)
            x=x_;
            j=j+1;
            if(z_<z_best)
                i
                z_best=z_
                x_best=x_;
                facLocs_best=facLocs_;
            end
         else
            if(acceptanceRule(delta,T,z_,z,thresholdParameter))
               x=x_;
               j=j+1;
            end
         end
         debug(end+1,:)=[z_best z z_ j L i];
      end
      tc=tc+changeTerminationCounter(j,L,fp);
      T=decreaseTemperature(j,L,ip,r,T);
      thresholdParameter=thresholdParameter*r;
      L=round(L*gamma)
      cycle=cycle+1
      cycles(cycle,:)=[z_best L T tc];
end
 
function result=thresholdAcceptance(delta,T,z_,z,mu)
result=z_<=(1+mu)*z;
 
function result=probabilisticAcceptance(delta,T,z_,z,mu)
result=exp(-delta/T)>rand;
 
function result=stoppingCondition(terminationCounter,cycles)
% result = terminationCounter >= 5; % alternative stopping condition
z_best=cycles(:,1);
if(length(z_best)<2)
    result=false;
    return;
end
change=(z_best(end)-z_best(end-1))/z_best(end);
eps_1=0.03;
if(change < eps_1)
    result=true;
else
    result=false;
end
 
function result=changeTerminationCounter(j,L,fp)
if(j/L <= fp)
   result=1;
else
   result=0;
end
 
function T=decreaseTemperature(j,L,ip,r,T)
if(j/L > ip)
   T=T/2;
else
   T=r*T;
end
 
function result=distance(x,y,degree)
if nargin < 3
    degree=2; % default: euclidean distance
end
result=(abs(x(1)-y(1))^degree+abs(x(2)-y(2))^degree)^(1/degree);
 
% neighborhood structure is based on the customer-facility assignments
% pick a neighbor of the current solution x where k is the number of changing assignments, p is the number of facilities
function x_=pickANeighbor(x,k,p)
changingCustomers=unidrnd(length(x),1,k);
x_=x;
for i=1:k
    newFacility=unidrnd(p-1);
    oldFacility=x(changingCustomers(i));
    if (newFacility >= oldFacility)
        newFacility=newFacility+1;
    end
    x_(changingCustomers(i))=newFacility;
end
 
function T=initializeTemperature(p_1,p,k,custLocs,demands)
n=100;
n_c=length(demands);
for i=1:n
    x=generateArbitrarySolution(n_c,p); % customer-facility assignments
    x_=pickANeighbor(x,k,p);
    facLocs=single_facility_optimization(x,custLocs,demands);
    facLocs_=single_facility_optimization(x_,custLocs,demands);
    Delta(i)=abs(f(x,demands,facLocs,custLocs)-f(x_,demands,facLocs_,custLocs));
end
T=mean(Delta)/log(1/p_1);
 
function x=generateArbitrarySolution(n_c,p)
x=unidrnd(p,1,n_c);
 
function result=findAssignments(facLocs,custLocs)
n_c=length(custLocs);
p=length(facLocs);
distances=distancesFromFacilities(facLocs,custLocs);
[minD,result]=min(distances');
 
% objective function value of a solution x
function result=f(x,demands,facLocs,custLocs)
n_c=length(custLocs);
result=0;
for cust=1:n_c
   facility=x(cust);
   dist=distance(facLocs(facility,:),custLocs(cust,:));
   result=result+dist*demands(cust);
end
 
function result=distancesFromSingleFacility(facilityLocation,customerLocations,degree)
if nargin < 3
    degree=2; % default: euclidean distance
end
for cust=1:length(customerLocations)
    y=customerLocations(cust,1:2);
    result(cust)=distance(facilityLocation,y,degree);
end
 
function distances=distancesFromFacilities(facLocs,custLocs)
n_c=length(custLocs);
for fac=1:length(facLocs)
    distances(1:n_c,fac)=distancesFromSingleFacility(facLocs(fac,:),custLocs);
end