Data Programming Course  Check-in [5cabac5e1f]

Many hyperlinks are disabled.
Use anonymous login to enable hyperlinks.

Overview
Comment:aggiunti gli script di soluzione delle vecchie lezioni e di dimostrazione usati
Downloads: Tarball | ZIP archive | SQL archive
Timelines: family | ancestors | descendants | both | trunk
Files: files | file ages | folders
SHA1:5cabac5e1f99b149e37edf0477504cb87e6372a5
User & Date: EnricoGiampieri 2017-03-16 15:47:50
Context
2017-03-17
09:17
terminate le slide dell'esercitazione finale check-in: 57209fe163 user: EnricoGiampieri tags: trunk
2017-03-16
15:47
aggiunti gli script di soluzione delle vecchie lezioni e di dimostrazione usati check-in: 5cabac5e1f user: EnricoGiampieri tags: trunk
15:37
inserito lo stub della seconda lezione extra su xarray check-in: 33ae3729a3 user: EnricoGiampieri tags: trunk
Changes

Added script_assortiti/MapReduce.m.

























































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
clc;
clear all;
close all;
% Big Data Analytics with Matlab

%==============  Datastore  =================
% Un datastore è un oggetto messo a disposizione
% da Matlab per la lettura di un singolo file o
% più dati. Il datastore agisce come una
% repository di dati che possiedono la stessa
% struttura e formato (ogni file nel datastore deve
% contenere dati dello stesso tipo e nello stesso
% "ordine", separati dallo stesso delimitatore.
%============================================

%============================================
% Creiamo ad esempio un datastore con i file di
% airline.csv, i quali includono partenze e arrivi
% di ogni compagnia aerea.
%============================================

ds = datastore('airlinesmall.csv');
ds

%============================================
% E' facile che nei database ci siano dei dati
% mancanti per cui è sempre meglio specificare
% come sono salvati, in modo da farli "trattare"
% in modo diverso.
%============================================

ds.TreatAsMissing = 'NA';

%============================================
% Da tutta la tabella dei dati assumiamo di
% concentrarci solamente su una variabile 
% e vediamo di fare un po' di statistica su
% questa.
%============================================

ds.SelectedVariableNames = {'Distance'};

%============================================
% Il bello dei datastore è che non è necessario
% caricarli tutti in memoria subito ma possiamo 
% dare anche solo una sbirciatina con i preview
%============================================

preview = preview(ds) % da notare l'assenza del ; 
                      % per farlo apparire in output

%============================================
% Ora possiamo lanciare il nostro MapReduce per il
% calcolo della media delle distanze tra voli.
%============================================

outds = mapreduce(ds, @MeanDistMapFun, @MeanDistReduceFun);

% e per vedere il risultato
readall(outds)

Added script_assortiti/MeanDistMapFun.m.























































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
%======================================
% La funzione Map riceve un chunk di dati,
% li riorganizza ed esegue eventuali calcoli
% preliminari prima di applicare la funzione
% per cui è effettivamente ideata.
% A livello pratico, in una sommatoria, procederà
% a sommare tra loro tutti le coppie chiave-valore
% in un dato intermedio detto KeyValueStore.
% Il numero di chiamate della funzione MAP da
% parte del MapReduce sarà pari al numero di chunks
% il cui si suddivide l'input.
% Il KeyValueStore intermedio raggruppa tutti i
% valori con chiavi uniche!
%======================================

function MeanDistMapFun( data, info, intermKVStore )
% Funzione MAP del paradigma Map Reduce per il calcolo
% della media di distanza tra i voli. 

% Inputs: data - dati effettivi
%         info - informazioni per la lettura legate al chunk
%         intermKVStore - KeyValueStore intermedio
%                         oggetto di cui la Map ha bisogno
%                         per sommare le coppie chiave-valore

% Considero solamente gli elementi che non sono vuoti
% all'interno dei miei dati
distances = data.Distance(~isnan(data.Distance));

% Creo un vettore di due elementi: distanza totale 
%                                  conteggio del chunk
sumLenValue = [sum(distances) length(distances)];

% Sommo i valori ottenuti nella variabile intermedia
% mettendogli la chiave 'sumAndLength'
add(intermKVStore, 'sumAndLength', sumLenValue);

% Lanciando il Map su tutti i chunk del data set
% avrò al termine la distanza totale e il conteggio
% per ogni chunk all'interno di intermKVStore

end

Added script_assortiti/MeanDistReduceFun.m.











































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
%======================================
% La funzione Reduce viene chiamata per
% ognuna delle chiavi uniche della Map.
% Ogni chiave può avere più valori salvati
% che vengono passati come ValueIterator.
% Per funzionare sono necessarie le funzioni
% 'hasnext' e 'getnext' per l'iterazione.
% Al termine vengono aggregati i valori intermedi
% e salvati in un'unica coppia chiave-valore.
%======================================

function MeanDistReduceFun( intermKey, intermValIter, outKVStore )
% Funzione di REDUCE per il paradigma Map Reduce
% per il calcolo della media delle distanze tra voli

% Inputs : intermKey - chiave intermedia (ogni chiamata specifia
%                       una nuova chiave unica)
%            intermValIter - interatore associato alla chiave attiva
%                           (contiene tutti i valori associati alla chiave attiva)
%                           Il passaggio tra i vari valori si esegue con la
%                           funzione 'hasnext'
%            outKVstore - KeyValueStore finale utilizzato per mettere
%                         insieme tutti i valori calcolati

% inizializzo a zero i valori
sumLen = [0 0];

% scorro tutte le coppie e sommo i valori intermedi
while hasnext(intermValIter)
    sumLen = sumLen + getnext(intermValIter);
end

% risommo tutto e metto tutto dentro la media
add(outKVStore, 'Mean', sumLen(1)/sumLen(2));

end

Added script_assortiti/PCA.m.

































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
clc;
clear all;
close all;

%% Load Dataset

load fisheriris;

gscatter(meas(:,1), meas(:,2), species,'rgb','osd');
xlabel('Sepal length');
ylabel('Sepal width');
N = size(meas,1);

data = meas;
label = species;
%% PCA

LblU = unique(label);
ColStr = {'k' 'r' 'g' 'b' 'm' 'c'};
tic;
Pca = princomp(data');
toc
Ipca = [1 2];
figure; hold on; box on
for i_l = 1:length(LblU)
    ii = strcmp(label, LblU(i_l));
    plot(Pca(ii,Ipca(1)), Pca(ii,Ipca(2)), ['o' ColStr{i_l}], 'linewidth', 2, 'markersize', 8)
end
legend(LblU, 'fontsize', 14);
axis square

%% FastPCA

data = data - mean(data);
h = 3;
sigma = cov(data);
col = length(data(1,:));
phi = rand(col,1);
p = 1;
eigen = zeros(col, h);
tic;
while(p < h)
   err = 1e4;
   while(err > 1e-8)
      phi = sigma * phi;
      if p > 1
          phi = phi - phi'*eigen(:,p-1);
      end
      phi = phi/norm(phi);
      eigen(:,p) = phi';
      err = abs(phi'*phi - 1);      
   end
   p = p + 1;
end
toc
fastPCA = data*eigen;

figure; hold on; box on
for i_l = 1:length(LblU)
    ii = strcmp(label, LblU(i_l));
    plot(fastPCA(ii,Ipca(1)), fastPCA(ii,Ipca(2)), ['o' ColStr{i_l}], 'linewidth', 2, 'markersize', 8)
end
legend(LblU, 'fontsize', 14);
axis square

Added script_assortiti/PCA.py.

















































































































































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 13 17:05:43 2016

@author: NICO
"""

from __future__ import division
import numpy as np
import matplotlib.pylab as plt
from sklearn import decomposition

"""
Funzione per la determinazione delle componenti principali (PCA) di una matrice di dati e label forniti in input. 
La funzione realizza il grafico delle prime due componenti principali e restituisce in output la matrice di queste.
"""

def PCA(Data, Label):
    Pca = decomposition.PCA()
    Pca.fit(Data)
    Pca.n_components = 2
    D_reduced = Pca.fit_transform(Data)
    
    return D_reduced
    

from sklearn import datasets


iris = datasets.load_iris()
data = iris.data
lbl = iris.target
Lbl = np.zeros(len(lbl))
lblU = np.unique(lbl)
for i in range(len(lbl)):
    if lbl[i] == lblU[0]:
        Lbl[i] = 0
    if lbl[i] == lblU[1]:
        Lbl[i] = 1
    if lbl[i] == lblU[2]:
        Lbl[i] = 2

plt.figure()
plt.subplot2grid((3,3+3), (0,1), colspan=2, rowspan=2)
plt.title('Original Data', fontsize=24, fontweight='bold')
plt.plot(data[:50,0], data[:50,1], 'bo', markersize=20)
plt.plot(data[100:,0], data[100:,1], 'rd', markersize=20)
plt.xlim(4,8.5)
plt.ylim(1, 5)
plt.subplot2grid((3,3+3), (0,0), rowspan=2)
plt.hist(data[:50,1], color='blue', orientation='horizontal', align='mid', alpha=0.5)
plt.hist(data[50:,1], color='red', orientation='horizontal', align='mid', alpha=0.5)
plt.ylabel("Feature 1", fontsize=20, fontweight='bold')
plt.ylim(1,5)
plt.subplot2grid((3,3+3), (2,1), colspan=2)
plt.hist(data[:50,0], color='blue', align='mid', alpha=0.5)
plt.hist(data[50:,0], color='red', align='mid', alpha=0.5)
plt.xlabel("Feature 2", fontsize=20, fontweight='bold')
plt.xlim(4, 8.5)

data2 = PCA(data, lbl)

#plt.figure()
plt.subplot2grid((3,3+3), (0,1+3), colspan=2, rowspan=2)
plt.title("Principal Component Analysis", fontsize=24, fontweight='bold')
plt.plot(data2[:50,0], data2[:50,1], 'bo', markersize=20)
plt.plot(data2[100:,0], data2[100:,1], 'rd', markersize=20)
plt.xlim()
plt.ylim()
plt.subplot2grid((3,3+3), (0,0+3), rowspan=2)
plt.hist(data2[:50,1], color='blue', orientation='horizontal', align='mid', alpha=0.5)
plt.hist(data2[50:,1], color='red', orientation='horizontal', align='mid', alpha=0.5)
plt.xlabel("2nd Component", fontsize=20, fontweight='bold')
plt.ylim()
plt.subplot2grid((3,3+3), (2,1+3), colspan=2)
plt.hist(data2[:50,0], color='blue', align='mid', alpha=0.5)
plt.hist(data2[50:,0], color='red', align='mid', alpha=0.5)
plt.xlabel("1st Component", fontsize=20, fontweight='bold')
plt.xlim()

#%%
#fastPCA algorithm

dd = data
dd -= np.mean(data, axis=0)
h = 2
sigma = np.cov(dd.T)
print sigma
d = len(data[0])
phi = np.random.rand(d)
print phi
p = 0
eigen = np.zeros((d, h))
while(p < h):
    err = 1e4
    while(err > 1e-16):
        phi = np.dot(sigma, phi)
        print "phi ", phi
        if p > 0:
            print "proj ", np.dot(phi.T, eigen[:,p-1])
            phi = phi - np.dot(phi.T, eigen[:,p-1])*eigen[:,p-1]
        phi /= np.linalg.norm(phi)
        print "phi gram ", phi
        eigen[:,p] = phi
        err = abs(np.dot(phi.T, phi)-1)
        print "err ", err
    p += 1

dd = np.dot(dd, eigen)

#%%
plt.figure()
plt.subplot(121)
plt.plot(data2[:50,0], data2[:50,1], 'bo', markersize=20)
plt.plot(data2[100:,0], data2[100:,1], 'rd', markersize=20)
plt.plot(data2[50:100,0], data2[50:100,1], 'gd', markersize=20)
plt.subplot(122)
plt.plot(dd[:50,0], dd[:50,1], 'bo', markersize=20)
plt.plot(dd[100:,0], dd[100:,1], 'rd', markersize=20)
plt.plot(dd[50:100,0], dd[50:100,1], 'gd', markersize=20)

Added script_assortiti/QDA.py.



















































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
# -*- coding: utf-8 -*-
"""
Created on Sat Mar 11 17:21:03 2017

@author: NICO
"""

from __future__ import division, print_function
import numpy as np

#%% A Naive version

def QDA_byhand(training, test, lbl):
    lblU = np.unique(lbl)
    classifier = np.zeros(len(lblU))
    for classe in range(0, len(lblU)):
        cls = training[lbl==lblU[classe],:]
        mu = np.mean(cls, 0)
        scatter = np.cov(cls.T)
        normalizer = 1/(2*np.pi*np.sqrt(np.linalg.det(scatter)))
        exponential = -0.5*np.dot((np.dot((test-mu),np.linalg.inv(scatter))), (test-mu).T)
        classifier[classe] = normalizer*np.exp(exponential)
        
    return lblU[classifier==max(classifier)]

#%% A more sofisticated version

def QDA_conditional(Train, Test, Label, conditional = 1e-6, prior = 1.):
    lblU = np.unique(Label)
    cls = np.empty(len(lblU))
    for i in range(len(lblU)):
        data = Train[Label == lblU[i], :].T
        diff =  Test - np.mean(data, axis = 1)
        cov = np.cov(data) + np.diag(np.full(len(data), conditional))
        #cls[i] = 1/(2*np.pi*np.sqrt(np.linalg.det(cov))) * np.exp(-0.5 *np.einsum('...k,kl,...l->...', Test-mean, np.linalg.inv(cov), Test-mean))
        cls[i] = -0.5 * np.einsum('...k,kl,...l->...', diff, np.linalg.inv(cov), diff) - 0.5 * np.log(np.linalg.det(cov)) + np.log(prior)#versione logaritmica ->meno costosa!
    
    return lblU[cls == cls.max()]


#%% Test of functions

from sklearn.qda import QDA
from sklearn.cross_validation import LeaveOneOut
from sklearn import datasets

iris = datasets.load_iris()
data = iris.data
label = iris.target

Loo = LeaveOneOut(len(label))
classifier = QDA()
Classes_SKlearn = np.empty(len(label))
Classes_byhand = np.empty(len(label))
Classes_conditional = np.empty(len(label))

for train, test in Loo:
    Training = data[train, :]
    Test = data[test, :]
    classifier.fit(Training, label[train])
    Classes_SKlearn[test] = classifier.predict(Test)
    Classes_byhand[test] = QDA_byhand(Training, Test, label[train])
    Classes_conditional[test] = QDA_conditional(Training, Test, label[train])

print ("Score QDA Sklearn %d over %d = %f%%"%(len(Classes_SKlearn[Classes_SKlearn == label]), 
                                      len(label), 
                                      float(len(Classes_SKlearn[Classes_SKlearn == label]))/len(label)))
print ("Score QDA by hand %d over %d = %f%%"%(len(Classes_byhand[Classes_byhand == label]), 
                                      len(label), 
                                      float(len(Classes_byhand[Classes_byhand == label]))/len(label)))
print ("Score QDA conditional %d over %d = %f%%"%(len(Classes_conditional[Classes_conditional == label]), 
                                      len(label), 
                                      float(len(Classes_conditional[Classes_conditional == label]))/len(label)))

Added script_assortiti/SparseMat.m.





































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
clc;
clear all;
close all;
% Sparse Matrix computation

% Iniziamo considerando una matrice FULL come la seguente

A = [ 0   0   0   5
      0   2   0   0
      1   3   0   0
      0   0   4   0]
% Al suo interno sono presenti numerosi valori nulli
% che sarebbe comodo non doversi portare dietro per
% velocizzare e snellire il proprio codice.

% Matlab mette a disposizione la funzione SPARSE
% che permette la conversione di una matrice 'densa'
% in una 'sparsa'

S = sparse(A)

% In alternativa, per creare direttamente una matrice
% sparsa è possibile fornire una coppia di vettori
% contenenti gli indici I e J, ed un vettore di Valori
% seguiti dalla dimensione N ed M della matrice
%
%               sparse(i, j, val, N, M)
%
% Nell'esempio di prima

S = sparse([3 2 3 4 1],[1 2 2 3 4],[1 2 3 4 5],4,4)


% Per tornare indietro alla definizione FULL basterà
% il comando

A = full(S)

%% Efficienza di calcolo
clc;
dim = 2*1e3;

% Creo una matrice random SPARSA con una percentuale di elementi non-nulli
A = sprand(dim, dim, 0.25);
tic, A*A; toc

Af = full(A);
tic, Af*Af; toc

%% Aumento il numero di zeri

A = sprand(dim, dim, 0.005);
Af = full(A);

tic, A*A; toc

tic, Af*Af; toc

%% Calcolo degli autovalori-autovettori

tic, [eigenval, eigenvec] = eigs(A); toc
tic, [eigenval, eigenvec] = eig(Af); toc




Added script_assortiti/metropolis.py.







































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
# -*- coding: utf-8 -*-
"""
Created on Sat Feb 11 18:37:12 2017

@author: NICO
"""
from __future__ import print_function
import numpy as np
import time

# Algorithm C++-like

N = 100000 # number of MC events
N_run = 100 # number of runs
Nhits = 0.0 # number of points accepted
pi = np.zeros(N_run) # values of pi

start_time = time.time() # start clock 
for I in range(N_run):
    Nhits = 0.0
    for i in range(N):
        x = np.random.rand()*2-1
        y = np.random.rand()*2-1
        res = x*x + y*y
        if res < 1:
            Nhits += 1.0
    pi[I] += 4. * Nhits/N

run_time = time.time()

print ("pi with ", N, " steps for ", N_run, " runs is ", np.mean(pi), " in ", run_time-start_time, " sec")



#%% Vectorized algorithm with Numpy functions

start_time = time.time() # start clock
pi = 4.*sum(np.where((np.random.rand(N_run, N)*2-1)**2+(np.random.rand(N_run, N)*2-1)**2 < 1, 1, 0))/N_run
run_time = time.time()
print ("pi with ", N, " steps for ", N_run, " runs is ", np.mean(pi), " in ", run_time-start_time, " sec")

#%% A more complex vectorization

start_time = time.time()
pi = np.mean((np.sum(np.random.uniform(-1, 1, [2, N_run, N])**2, axis=0)<1), axis=1)*4
run_time = time.time()
print ("pi with ", N, " steps for ", N_run, " runs is ", np.mean(pi), " in ", run_time-start_time, " sec")