%
%   Copyright (C) %{YEAR} by %{AUTHOR}
%   <%{EMAIL}>
%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%

% Ritorna i Bit Plabes dell'immagine a toni di grigio

function showBitPlanes(img)

    imgGray = double( rgb2gray(img) );
    titleString = 'bit planes ';
    
    % MSB ... LSB
    k = 128;
    
    for b=1:8
        
        subplot(2, 4, b);
        imshow( (bitand(imgGray, k) / k) * 255 ); % Fa un and dei bit
        title([titleString int2str(b-1)]);
        k = k/2; % Shifta di 2 i bit
        
    end;
    
return;

--- 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

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

function [ salida ] = pasabajos(nombre)
%Filtro pasabajos de 1er orden:

[X,Fs,bits] = wavread(nombre); %carga el wav
FX = fft(X);
Xmed=mean(X); %debe dar 0
Xvar=var(X);

pfiltro = [1 1]; %orden 1, promediador de 2 muestras: y[n]=(x[n]+x[n-1])/2

Y=1/2*filter(pfiltro,[1],X);
FY = fft(Y);

N=length(FX);
if ( mod(N,2) == 0 )
t=N/2;
else
t=floor(N/2)+1;
end
rectaFrec=linspace(0,pi*Fs/(2*pi),floor(N/2)+1);

FX = fftshift(FX);
FY = fftshift(FY);

subplot(2,1,1)
title(’Original’)
Espectro=abs(FX(t:N));
plot(rectaFrec,Espectro,’r'); %Grafico en funcion de la frecuencia en Hz

subplot(2,1,2)
title(’Pasabajos de 1er orden’)
Espectro=abs(FY(t:N));
plot(rectaFrec,Espectro,’b'); %Grafico en funcion de la frecuencia en Hz

salida = Y;

--- do_experiments.m
% Usage of heuristics in console:
% [p,l,d]=shortest_path('d198.txt');
% [p2,l2,d2]=insertion_heuristics('d198.txt',@farthest_insertion);
% [p2,l2,d2]=insertion_heuristics('d198.txt',@arbitrary_insertion);

function results=do_experiments()
data_files={'d198.txt';'d1291.txt';'rat575.txt'};
number_of_experiments=10;
results=[];
for idx_data=1:length(data_files)
    for i=1:number_of_experiments
        data=char(data_files(idx_data));
        [path,total_length(1),dist]=shortest_path(data);
        [path,total_length(2)]=opt2(path,dist);
        [path,total_length(3),dist]=insertion_heuristics(data,@farthest_insertion);
        [path,total_length(4)]=opt2(path,dist);
        [path,total_length(5),dist]=insertion_heuristics(data,@arbitrary_insertion);
        [path,total_length(6)]=opt2(path,dist);
        results(i+(idx_data-1)*number_of_experiments,:)=total_length;
    end
end

--- shortest_path.m
% calculates shortest path using nearest neighborhood
function [path,total_length,distances] = shortest_path(data_file,idx_initial_city)
distances=distance(read_location_data(data_file));
n=length(distances);
if nargin < 2
    idx_initial_city=ceil(rand(1)*n);
end
path=[];
outside=(1:n);
[path,outside]=add_to_path(path,outside,outside(idx_initial_city));
for i = 1:n-1
    entering_city=nearest_neighbor(distances,path(end),outside);
    [path,outside]=add_to_path(path,outside,entering_city);
end
total_length=total_length_of_cycle(distances,path);

function [path,outside] = add_to_path(path,outside,city)
path(end+1)=city;
outside(outside==city)=[];

--- distance.m
function d = distance(a)
n=length(a);
d=zeros(n,n);
for i = 1:n
    for j = i+1:n
        d(i,j)=norm(a(i,:)-a(j,:));
        d(j,i)=d(i,j);
    end;
end;

--- read_location_data.m
function positions = read_location_data(data_file)
load(data_file);
filename=regexp(data_file,'\w*(?=\.)','match');
positions=eval(filename{1});

--- nearest_neighbor.m
function nearest_city = nearest_neighbor(distances,last_city,outside)
    min_distance = inf;
    for j=1:length(outside)
        neighbors_distance = distances(outside(j),last_city);
        if neighbors_distance < min_distance
            min_distance = neighbors_distance;
            nearest_city = outside(j);
        end
    end
end
    
--- total_length_of_cycle.m
function total_length = total_length_of_cycle(distances,path)
total_length = 0;
for i = 1:length(path)-1
    total_length = total_length + distances(path(i),path(i+1));
end
total_length = total_length + distances(path(end),path(1));

--- insertion_heuristics.m
% generic insertion algorithm
function [path,total_length,distances] = insertion_heuristics(data_file,insertion_rule,idx_initial_city)
distances=distance(read_location_data(data_file));
n=length(distances);
if nargin < 3
    idx_initial_city=ceil(rand(1)*n);
end
path=[];
outside=(1:n);
[path,outside]=add_to_path(path,outside,outside(idx_initial_city),1);
for i = 1:n-1
    entering_city=insertion_rule(path,outside,distances);
    entry_position=find_entry_position(distances,entering_city,path);
    [path,outside]=add_to_path(path,outside,entering_city,entry_position);
end
total_length=total_length_of_cycle(distances,path);

function [path,outside] = add_to_path(path,outside,entering_city,entry_position)
old_path=path;
outside(outside==entering_city)=[];
path(entry_position)=entering_city;
for i=entry_position:length(old_path)
    path(i+1)=old_path(i);
end

function entry_position=find_entry_position(distances,entering_city,path)
min_increase_in_length = inf;
global path_length
path_length = length(path);
for i=1:path_length
    before = distances(path(i),path(r(i+1)));
    after = distances(entering_city,path(i))+distances(entering_city,path(r(i+1)));
    increase_in_length = after-before;
    if increase_in_length < min_increase_in_length
        min_increase_in_length = increase_in_length;
        entry_position = r(i+1);
    end
end
    
% if index is greater than the path length, turn the index for one round
function result = r(index)
global path_length
if index > path_length
    result = index - path_length;
else
    result = index;
end

--- farthest_insertion.m
function entering_city=farthest_insertion(path,outside,distances)
max_distance=0;
for i=1:length(path)
    for j=1:length(outside)
        if distances(path(i),outside(j)) > max_distance
            max_distance=distances(path(i),outside(j));
            entering_city=outside(j);
        end
    end
end

--- arbitrary_insertion.m
function entering_city = arbitrary_insertion(path,outside,distances)
entering_city=outside(ceil(rand(1)*length(outside)));

--- opt2.m
function [path,total_length] = opt2(path,distances)
global n 
n = length(path);
for i=1:n
    for j=1:n-3
        if change_in_path_length(path,i,j,distances) < 0 
            path=change_path(path,i,j);
        end
    end
end
total_length = total_length_of_cycle(distances,path);
        
function result = change_in_path_length(path,i,j,distances)
before=distances(path(r(i)),path(r(i+1)))+distances(path(r(i+1+j)),path(r(i+2+j)));
after=distances(path(r(i)),path(r(i+1+j)))+distances(path(r(i+1)),path(r(i+2+j)));
result=after-before;
            
function path = change_path(path,i,j)
old_path=path;
% exchange edge cities
path(r(i+1))=old_path(r(i+1+j));
path(r(i+1+j))=old_path(r(i+1));
% change direction of intermediate path 
for k=1:j-1
    path(r(i+1+k))=old_path(r(i+1+j-k));
end

% if index is greater than the path length, turn index one round off
function result = r(index)
global n
if index > n
    result = index - n;
else
    result = index;
end