% Questo programma consiste di due parti, la prima riordina un set di n numeri usando un metodo
% basato sul confronto tra coppie (ordine quadratico) ed un secondo metodo che sfrutta
% il partizionamento (ordine nlog(n)). Questo viene fatto per diversi valori di n.
% Nella seconda parte, per ciascun metodo vengono stimati i parametri di un modello teorico
% per il costo computazionale del metodo in funzione di n e vengono plottati i due modelli
% assieme ai costi misurati.
% Commenti a cura di Luigi Lunardon e Vittorio Maio
clear
close all
m=10;
for jl=1:m
n=100*jl
num(jl)=n;
x=rand(1,n);
zx=x;
% In questo specifico caso riordineremo m vettori creati in maniera random formati da un numero di elementi che
% varia da 100 a 100*m. Nel vettore zx salviamo il vettore x in modo da utilizzare per ogni n
% lo stesso set di numeri con entrambi i metodi.
% Metodo basato sul confronto tra coppie
tic % Calcoliamo il costo computazionale
for i=1:n-1
imin1=find(x==min(x));
imin=min(imin1);
y(i)=min(x);
w=ones(length(x)-1,1);
kk=1;
for k=1:length(x)
if k == imin
else
w(kk)=x(k);
kk=kk+1;
end
end
x=w;
end
y(n)=x;
tempo1(jl)=toc; % salviamo in tempo1 il tempo di calcolo
x=zx;
% Metodo basato sul partizionamento
tic % Calcoliamo il costo computazionale
% L'idea del metodo di partizionamento e' quella di operare ricorsivamente la divisione
% dei numeri di un insieme in due sottinsiemi: da un lato gli elementi <= di un certo
% valore e nell'altro quelli > (in questo programma l'elemento di confronto e' la media
% aritmetica dei numeri). Il costo computazionale medio e' ridotto a n*log(n). Il codice seguente
% utilizza la funzione part, da noi preparata e salvata col nome di part.m. L'input di part
% e' un vettore x di numeri. L'output e' un vettore xx con le stesse dimensioni, le cui
% prime k componenti contengono le componenti di x minori o uguali alla soglia, mentre le
% restanti componenti sono le rimanenti componenti di x. In output part fornisce anche il valore
% di k. Data la partizione corrente degli n numeri, quella successiva viene costruita applicando
% la funzione part a ciascun elemento della partizione.
[xx nn]=part(x);
k=2; % numero iniziale di elementi della partizione
kold=1;
index=zeros(1,3);
index(2)=nn;
index(3)=n; % posizioni degli estremi sinistri-1 degli elementi della partizione corrente.
% index e' un vettore la cui componente i indica che l'elemento i della partizione corrente
% ha inizio nella posizione i+1. L'ultima componente di index e' sempre n. Il metodo ha bisogno
% di un criterio di arresto. Percio' non possiamo usare un ciclo for, ma usiamo un ciclo while.
% Il ciclo while in matlab lavora esattamente come in C, all'inizio del ciclo verifica una certa
% condizione e se e' vera continua il ciclo. Il criterio d'arresto e' qui che la partizione prima e
% dopo un'iterazione non cambi numero di elementi.
while k>kold
kold=k;
index1=zeros(1);
x=xx;
kk=1; % contatore degli elementi della partizione successiva
for j=1:k % loop sugli elementi della partizione corrente
z=x(index(j)+1:index(j+1)); % vettore dell'elemento j della partizione corrente
if max(z)>min(z) %Se z non ha gli elementi tutti uguali, viene applicata part
[zz nn]=part(z);
index1(kk+1)=index(j)+nn;
index1(kk+2)=index(j+1); % posizioni degli estremi sinistri-1 degli elementi della partizione successiva
kk=kk+2; % abbiamo diviso l'elemento j della partizione corrente in due
else
zz=z;
index1(kk+1)=index(j+1); % non abbiamo diviso l'elemento j della partizione corrente in due
kk=kk+1;
end
xx(index(j)+1:index(j+1))=zz; % Sostituiamo il vettore z dell'elemento j della partizione corrente con quello diviso in due (zz).
end
k=kk-1; % nuovo numero di elementi della partizione
index=index1;
end
tempo2(jl)=toc; % salviamo in tempo2 il tempo di calcolo
cc(jl)= sum(abs(xx-y));
end
sum(cc) % verifica che i due metodi danno lo stesso risultato
plot(num,tempo1,'+')
hold on
plot(num,tempo2,'*')
% Abbiamo plottato sullo stesso grafico i tempi corrispondenti al riordinamento di un vettore
% prima con il primo metodo e poi con il secondo. Si tratta ora di determinare il modello teorico
% del tempo impiegato versus n. Per fare cio' useremo un fit non lineare ossia un'interpolazione dei dati
% sperimentali che sfrutta la conoscenza a priori del modello teorico non lineare (in questo caso
% in realta' i due modelli sono lineari, ma utilizziamo lo stesso il fit non lineare per scopi
% didattici). Vediamo come funziona il comando nlinfit: nlinfit(indip, dip, @model, beta).
% indip e dip sono rispettivamente i vettori con i valori della variabile indipendente e
% dipendente della funzione non lineare. model e' una funzione matlab da noi definita con cui, %a partire dal vettore dei valori della variabile indipendente e di quello dei parametri (beta),
% si ottiene il vettore dei valori della variabile dipendente. Il metodo di minimizzazione
% alla base di nlinfit e' locale e richiede l'inizializzazione dei valori dei parametri (vettore
% beta). La funzione nlinfit da' come output il vettore beta con i suoi valori ottimali (quelli
% che minimizzano la somma degli scarti tra i valori sperimentali della variabile dipendente
% e quelli calcolati col modello teorico a partire dai valori della variabile indipendente).
beta(1)=tempo2(m)/num(m)/log(num(m)); % inizializzazione del parametro del modello
beta=nlinfit(num, tempo2, @model2, beta);
zzz=linspace(num(1), num(m), 1000);
yyy=model2(beta, zzz);
beta(1)=tempo1(m)/num(m)^2;
beta(2)=0; % inizializzazione dei parametri del modello
beta=nlinfit(num, tempo1, @model1, beta);
xxx=model1(beta, zzz);
% Calcolate le interpolanti non resta che plottarle.
plot(zzz, xxx)
plot(zzz, yyy)
print -dpdf plot1.pdf
%--------------------------------------------------------------------------------------
% mettere in un file chiamato part.m
% Quella di seguito e' una funzione MATLAB da noi costruita. Le funzioni sono sequenze di
% comandi esterne al main e strutturate per adempiere ad un compito specifico. Le funzioni vengono
% scritte, salvate come file nome_funzione.m e poi richiamate all'interno del main, scrivendo
% output=nome_funzione(input). Le variabili interne alla funzione sono chiamate variabili
% locali e non sono visibili nel main se non sono di output. La funzione qui di seguito ha in
% input un vettore x, ed in output un vettore y ottenuto come partizione di x.
% Le prime k componenti di y sono i valori di x minori o uguali della media delle componenti di x
% mentre le restanti componenti hanno valore maggiore della media. In output abbiamo anche
% il valore di k.
% La sintassi usata per definire una funzione e' la seguente: %function [valori di output non separati da virgole]=nome_funzione(elementi dell'imput separati
% da virgole)
function [y n1] = part(x)
m=min(x);
M=max(x);
% Le funzioni max(x) e min(x) restituiscono rispettivamente il massimo e il minimo valore di un
% vettore x.
p=(M+m)/2;
if m < M % non facciamo nessuna operazione se gli elementi di x sono tutti uguali fra loro.
% Altrimenti per prima cosa salviamo nel vettore i1 le posizioni delle componenti di x che sono
% minori o uguali a p. Per questo utilizziamo il comando find(x<=p). Il comando find ha come
% argomento la condizione richiesta e crea il vettore con le posizioni delle componenti
% che la soddisfano.
i1=find(x<=p);
x1=x(i1); % Scriviamo il vettore x1 che ha come elementi gli elementi di x che
% si trovano nelle posizioni indicate da i1.
n1=length(x1);
i2=find(x>p);% Ripetiamo lo stesso procedimento per creare il vettore x2
x2=x(i2);
y=[ x1 x2]; % scriviamo il vettore y come concatenazione dei due.
else
y=x;
'ATTENZIONE dati tutti uguali'
end
%--------------------------------------------------------------------------------------
% mettere in un file chiamato model1.m
function yy=model(beta,ttt)
yy=beta(1)*ttt.^2+beta(2)*ttt;
%--------------------------------------------------------------------------------------
% mettere in un file chiamato model2.m
function yy=model(beta,ttt)
yy=beta(1)*ttt.*log(ttt);