mp4: done 1a-c in MATLAB

mp4/Makefile

@@ -0,0 +1,15 @@
+filename=template
+
+pdf:
+	pdflatex ${filename}
+	#bibtex ${filename}
+	pdflatex ${filename}
+	pdflatex ${filename}
+	make clean
+
+read:
+	evince ${filename}.pdf &
+
+
+clean:
+	rm -f *.out *.log *.bbl *.blg *.aux ${filename}.log  ${filename}.ps ${filename}.aux ${filename}.out ${filename}.dvi ${filename}.bbl ${filename}.blg

mp4/Project_4_Maggioni_Claudio.tex

@@ -0,0 +1,30 @@
+\documentclass[unicode,11pt,a4paper,oneside,numbers=endperiod,openany]{scrartcl}
+
+\input{assignment.sty}
+
+
+
+
+\begin{document}
+
+
+\setassignment
+\setduedate{Wednesday, 18 November 2020, 11:55 PM}
+
+\serieheader{Numerical Computing}{2020}{Student: FULL NAME}{Discussed with: FULL NAME}{Solution for Project 4}{}
+\newline
+
+\assignmentpolicy
+
+
+\begin{enumerate}
+
+\item  \textbf{Spectral clustering of non-convex sets [60 points]:}
+
+
+\item \textbf{Spectral clustering of real-world graphs [40 points]:}
+
+\end{enumerate}
+
+
+\end{document}

mp4/Project_4_Maggioni_Claudio/datasets/clusterincluster.m

@@ -0,0 +1,43 @@
+function data = clusterincluster(N, r1, r2, w1, w2, arms)
+
+    if nargin < 1
+        N = 1000;
+    end
+    if nargin < 2
+        r1 = 1;
+    end
+    if nargin < 3
+        r2 = 5*r1;
+    end
+    if nargin < 4
+        w1 = 0.8;
+    end
+    if nargin < 5
+        w2 = 1/3;
+    end
+    if nargin < 6
+        arms = 64;
+    end
+    
+    data = [];
+    
+    N1 = floor(N/2);
+    N2 = N-N1;
+    
+    phi1 = rand(N1,1) * 2 * pi;
+%   dist1 = r1 + randint(N1,1,3)/3 * r1 * w1;
+    dist1 = r1 + (randi(3,N1,1)-1)/3 * r1 * w1;
+    d1 = [dist1 .* cos(phi1) dist1 .* sin(phi1) zeros(N1,1)];
+
+    perarm = round(N2/arms);
+    N2 = perarm * arms;
+    radperarm = (2*pi)/arms;
+    phi2 = ((1:N2) - mod(1:N2, perarm))/perarm * (radperarm);
+    phi2 = phi2';
+    dist2 = r2 * (1 - w2/2) + r2 * w2 * mod(1:N2, perarm)'/perarm;
+    d2 = [dist2 .* cos(phi2) dist2 .* sin(phi2) ones(N2,1)];    
+    
+    data = [d1;d2];   
+
+    scatter(data(:,1), data(:,2), 20, data(:,3)); axis square;
+end

mp4/Project_4_Maggioni_Claudio/datasets/corners.m

@@ -0,0 +1,35 @@
+function data = corners(N, scale, gapwidth, cornerwidth)
+
+    if nargin < 1
+        N = 1000;
+    end
+    if mod(N,8) ~= 0
+        N = round(N/8) * 8;
+    end
+
+    if nargin < 2
+        scale = 10;
+    end
+    if nargin < 3
+        gapwidth = 2;
+    end   
+    if nargin < 4
+        cornerwidth = 2;
+    end
+
+    perCorner = N/4;
+
+    xplusmin = [ones(perCorner,1); -1*ones(perCorner,1); ones(perCorner,1); -1*ones(perCorner,1)];
+    yplusmin = [ones(perCorner,1); -1*ones(2*perCorner,1); ones(perCorner,1)];
+    
+    horizontal = [xplusmin(1:2:end) * gapwidth + xplusmin(1:2:end) * scale .* rand(N/2,1), ...
+                  yplusmin(1:2:end) * gapwidth + cornerwidth * yplusmin(1:2:end) .* rand(N/2,1), ...
+                  floor((0:N/2-1)'/(perCorner*.5))];
+       
+    vertical = [xplusmin(2:2:end) * gapwidth + cornerwidth * xplusmin(2:2:end) .* rand(N/2,1), ...
+                yplusmin(2:2:end) * gapwidth + yplusmin(2:2:end) * scale .* rand(N/2,1), ...
+                floor((0:N/2-1)'/(perCorner*.5))];
+    
+    data=  [horizontal/2; vertical/2];
+
+end

mp4/Project_4_Maggioni_Claudio/datasets/crescentfullmoon.m

@@ -0,0 +1,31 @@
+function data = crescentfullmoon(N, r1, r2, r3)
+
+if nargin < 1
+    N = 1000;
+end
+if mod(N,4) ~= 0
+    N = round(N/4) * 4;
+end
+if nargin < 2
+    r1 = 5;
+end
+if nargin < 3
+    r2 = 10;
+end
+if nargin < 4
+    r3 = 15;
+end
+
+N1 = N/4;
+N2 = N-N1;
+
+phi1 = rand(N1,1) * 2 * pi;
+R1 = sqrt(rand(N1, 1));
+moon = [cos(phi1) .* R1 * r1 sin(phi1) .* R1 * r1 zeros(N1,1)];
+
+d = r3 - r2;
+phi2 = pi + rand(N2,1) * pi;
+R2 = sqrt(rand(N2,1));
+crescent = [cos(phi2) .* (r2 + R2 * d) sin(phi2) .* (r2 + R2 * d) ones(N2,1)];
+
+data = [moon/2; crescent/2];

mp4/Project_4_Maggioni_Claudio/datasets/getGraphs.m

@@ -0,0 +1,30 @@
+function [Barth,Airfoil,Crack] = getGraphs()
+% function to get the Adjacency matrix and the coordinate list of various graphs
+% dimosthenis.pasadakis@usi.ch
+% ICS, USI.
+clear all;
+close all;
+
+
+addpath Meshes
+% 
+load('barth.mat');
+Barth.W   = Problem.A;
+Barth.Pts = Problem.aux.coord;
+% 
+load('airfoil1.mat');
+Airfoil.W   = Problem.A;
+Airfoil.Pts = Problem.aux.coord;
+% 
+load('crack_dual.mat');
+Crack.W   = Problem.A;
+Crack.Pts = Problem.aux.coord;
+
+figure;
+subplot(131)
+gplotLength(Barth.W,Barth.Pts,2);
+subplot(132)
+gplotLength(Airfoil.W,Airfoil.Pts,2);
+subplot(133)
+gplotLength(Crack.W,Crack.Pts,2);
+end

mp4/Project_4_Maggioni_Claudio/datasets/getPoints.m

@@ -0,0 +1,63 @@
+function [Pts_spirals,Pts_clusterin,Pts_corn,Pts_halfk,Pts_fullmoon,Pts_out] = getPoints()
+% Function to get the coordinate list of various pointclouds 
+% ----------------------------------------------------------
+% Copyright (c) 2013, Jeroen Kools
+% 
+% Modified for academic purposes by D. Pasadakis
+% dimosthenis.pasadakis@usi.ch
+% ICS, USI.
+
+clear all;
+close all;
+
+figure;
+hold on;
+dotsize = 12;
+ colormap([1 0 .5;   % magenta
+           0 0 .8;   % blue
+           0 .6 0;   % dark green
+           .3 1 0]); % bright green
+
+subplot(231);
+Pts_spirals = twospirals();
+scatter(Pts_spirals(:,1), Pts_spirals(:,2), dotsize,'k'); axis equal;
+% scatter(Pts_spirals(:,1), Pts_spirals(:,2), dotsize, Pts_spirals(:,3)); axis equal; % visualization 
+% axis off;
+title('Two spirals');
+
+subplot(232);
+Pts_clusterin = clusterincluster();
+scatter(Pts_clusterin(:,1), Pts_clusterin(:,2), dotsize,'k'); axis equal;
+% scatter(Pts_clusterin(:,1),Pts_clusterin(:,2), dotsize, Pts_clusterin(:,3)); axis equal; % visualization 
+% axis off;
+title('Cluster in cluster');
+
+subplot(233);
+Pts_corn = corners();
+scatter(Pts_corn(:,1), Pts_corn(:,2), dotsize,'k'); axis equal;
+% scatter(Pts_corn(:,1),Pts_corn(:,2), dotsize, Pts_corn(:,3)); axis equal; % visualization 
+% axis off;
+title('Corners');
+
+subplot(234);
+Pts_halfk = halfkernel();
+scatter(Pts_halfk(:,1), Pts_halfk(:,2), dotsize,'k'); axis equal;
+% scatter(Pts_halfk(:,1),Pts_halfk(:,2), dotsize, Pts_halfk(:,3)); axis equal; % visualization 
+% axis off;
+title('Half-kernel');
+
+subplot(235);
+Pts_fullmoon = crescentfullmoon();
+scatter(Pts_fullmoon(:,1), Pts_fullmoon(:,2), dotsize,'k'); axis equal;
+% scatter(Pts_fullmoon(:,1),Pts_fullmoon(:,2), dotsize,Pts_fullmoon(:,3)); axis equal; % visualization 
+% axis off;
+title('Crescent & Full Moon');
+
+subplot(236);
+Pts_out = outlier();
+scatter(Pts_out(:,1), Pts_out(:,2), dotsize,'k'); axis equal;
+% scatter(Pts_out(:,1),Pts_out(:,2),dotsize,Pts_out(:,3)); axis equal; % visualization 
+% axis off;
+title('Outlier');
+
+end

mp4/Project_4_Maggioni_Claudio/datasets/gplotLength.m

@@ -0,0 +1,48 @@
+function [d] = gplotLength(a, x, linewidth)
+% function [d] = gplotFancy2(a, x,linewidth)
+%
+% plot with colors depending on lengths of
+% edges in x,
+% a: Adjacency
+% x: coords
+% linewidth: for ploting
+
+[ai,aj] = find(a);
+m = length(ai);
+d = double(a);
+
+for i = 1:m, 
+    d(ai(i),aj(i)) = norm(x(ai(i),:)-x(aj(i),:)); 
+end
+
+[di,dj,dv] = find(d);
+
+sdv = sort(dv);
+l = length(sdv);
+
+ran = [64:-1:0]/64;
+ran = ran.^2;
+ran = 1-ran;
+
+ll = round([1:l/65:l]);
+ranges = sdv(ll);
+
+hold on;
+set(gca,'Color',[0 0 0])
+
+cm = jet;
+
+%cm = .3+.7*cm;
+
+for i = 1:64;
+ind = (dv <= ranges(i+1))&(dv>ranges(i));
+
+thisa = sparse(ai(ind),aj(ind),1);
+[gpx,gpy] = gplot(thisa,x);
+if (~isempty(gpx)),
+  ph(i) = plot(gpx,gpy,'-','Color',cm(i,:),'LineWidth',linewidth);
+  axis off
+end
+
+end
+

mp4/Project_4_Maggioni_Claudio/datasets/halfkernel.m

@@ -0,0 +1,32 @@
+function data = halfkernel(N, minx, r1, r2, noise, ratio)
+
+if nargin < 1
+    N = 1000;
+end
+if mod(N,2) ~= 0
+    N = N + 1;
+end
+if nargin < 2
+    minx = -20;
+end
+if nargin < 3
+    r1 = 20;
+end
+if nargin < 4
+    r2 = 35;
+end
+if nargin < 5
+    noise = 4;
+end
+if nargin < 6
+    ratio = 0.6;
+end
+
+phi1 = rand(N/2,1) * pi;
+inner = [minx + r1 * sin(phi1) - .5 * noise  + noise * rand(N/2,1) r1 * ratio * cos(phi1) - .5 * noise + noise * rand(N/2,1) ones(N/2,1)];
+    
+phi2 = rand(N/2,1) * pi;
+outer = [minx + r2 * sin(phi2) - .5 * noise  + noise * rand(N/2,1) r2 * ratio * cos(phi2) - .5 * noise  + noise * rand(N/2,1) zeros(N/2,1)];
+
+data = [inner/4; outer/4];
+end

mp4/Project_4_Maggioni_Claudio/datasets/outlier.m

@@ -0,0 +1,38 @@
+function data = outlier(N, r, dist, outliers, noise)
+
+    if nargin < 1
+        N = 600;
+    end
+    if nargin < 2
+        r = 20;
+    end
+    if nargin < 3
+        dist = 30;
+    end
+    if nargin < 4
+        outliers = 0.04;
+    end
+    if nargin < 5
+        noise = 5;
+    end
+
+    N1 = round(N * (.5-outliers));
+    N2 = N1;
+    N3 = round(N * outliers);
+    N4 = N-N1-N2-N3;
+
+    phi1 = rand(N1,1) * pi;
+    r1 = sqrt(rand(N1,1))*r;
+    P1 = [-dist + r1.*sin(phi1) r1.*cos(phi1) zeros(N1,1)];
+
+    phi2 = rand(N2,1) * pi;
+    r2 = sqrt(rand(N2,1))*r;
+    P2 = [dist - r2.*sin(phi2) r2.*cos(phi2) 3*ones(N2,1)];    
+    
+    P3 = [rand(N3,1)*noise dist+rand(N3,1)*noise 2*ones(N3,1)];    
+    
+    P4 = [rand(N4,1)*noise -dist+rand(N4,1)*noise ones(N4,1)];
+    
+    data = [P1/4.5; P2/4.5; P3/4.5; P4/4.5];
+
+end

mp4/Project_4_Maggioni_Claudio/datasets/twospirals.m

@@ -0,0 +1,34 @@
+function data = twospirals(N, degrees, start, noise)
+% Generate "two spirals" dataset with N instances.
+% degrees controls the length of the spirals
+% start determines how far from the origin the spirals start, in degrees
+% noise displaces the instances from the spiral. 
+%  0 is no noise, at 1 the spirals will start overlapping
+
+    if nargin < 1
+        N = 2000;
+    end
+    if nargin < 2
+        degrees = 570;
+    end
+    if nargin < 3
+        start = 90;
+    end
+    if nargin < 5
+        noise = 0.2;
+    end  
+    
+    deg2rad = (2*pi)/360;
+    start = start * deg2rad;
+
+    N1 = floor(N/2);
+    N2 = N-N1;
+    
+    n = start + sqrt(rand(N1,1)) * degrees * deg2rad;   
+    d1 = [-cos(n).*n + rand(N1,1)*noise sin(n).*n+rand(N1,1)*noise zeros(N1,1)];
+    
+    n = start + sqrt(rand(N1,1)) * degrees * deg2rad;      
+    d2 = [cos(n).*n+rand(N2,1)*noise -sin(n).*n+rand(N2,1)*noise ones(N2,1)];
+    
+    data = [d1/2;d2/2];
+end

mp4/Project_4_Maggioni_Claudio/src/ClusterGraphs.m

@@ -0,0 +1,65 @@
+% Cluster 2D real-world graphs with spectral clustering and compare with k-means
+% USI, ICS, Lugano
+% Numerical Computing 
+
+clear all;close all;
+warning OFF;
+
+addpath ../datasets
+addpath ../datasets/Meshes
+
+load airfoil1.mat
+% load barth.mat
+% load grid2.mat
+% load 3elt.mat
+
+% Specify the number of clusters
+K = 4;
+% Read graph
+W   = Problem.A;
+Pts = Problem.aux.coord;
+n   = size(Pts,1);
+% dummy var
+dummy_map = ones(n,1);
+
+figure;
+spy(W)
+title('Adjacency matrix')
+%% 2a) Create the Laplacian matrix and plot the graph using the 2nd and 3rd eigenvectors
+% \----------------------------/
+%     Your implementation
+% \----------------------------/
+
+% Eigen-decomposition
+% \----------------------------/
+%     Your implementation
+% \----------------------------/
+
+% Plot and compare
+figure;
+subplot(1,2,1);
+gplot(W,Pts)
+xlabel('Nodal coordinates')
+subplot(1,2,2);
+gplot(W,Pts)
+xlabel('TODO: Plot using Eigenvector coordinates')
+
+%% 2b) Cluster each graph in K = 4 clusters with the spectral and the 
+%      k-means method, and compare yourresults visually for each case.
+
+% \----------------------------/
+%     Your implementation
+% \----------------------------/
+
+
+% Compare and visualize
+figure;
+subplot(1,2,1);
+gplotmap(W,Pts,dummy_map)
+title('TODO: Plot the spectral clusters')
+subplot(1,2,2);
+gplotmap(W,Pts,dummy_map)
+title('TODO: Plot the K-means clusters')
+
+%% 2c) Calculate the number of nodes per cluster
+[Spec_nodes,Kmeans_nodes] = USI_ClusterMetrics(K,dummy_map,dummy_map);

mp4/Project_4_Maggioni_Claudio/src/ClusterMetrics.m

@@ -0,0 +1,13 @@
+function [Spec_nodes,Kmeans_nodes] = ClusterMetrics(K,x_spec,x_kmeans)
+%% METRICS
+Spec_nodes   = 0;
+Kmeans_nodes = 0;
+
+ % 2c) Calculate the number of nodes per cluster (for k-means and spectral
+ %     clustering) and plot them in a histogram.
+ 
+% \----------------------------/
+%     Your implementation
+% \----------------------------/
+ 
+end

mp4/Project_4_Maggioni_Claudio/src/ClusterPoints.m

@@ -0,0 +1,64 @@
+% Cluster 2D pointclouds with spectral clustering and compare with k-means
+% USI, ICS, Lugano
+% Numerical Computing 
+
+clear variables;
+close all;
+warning OFF;
+
+addpath ../datasets
+addpath ../datasets/Meshes
+
+% Specify the number of clusters
+K = 2;
+
+%% 1a) Get coordinate list from pointclouds
+% Coords used in this script
+Pts = getPoints();
+figure;
+scatter(Pts(:,1),Pts(:,2))
+title('Two Spirals')
+
+n = size(Pts, 1);
+
+% Create Gaussian similarity function
+[S] = similarityfunc(Pts(:,1:2));
+
+%% 1b) Find the minimal spanning tree of the full graph. Use the information
+%      to determine a valid value for epsilon        
+H = minSpanTree(S);
+epsilon = max(H(H > 0), [], 'all');
+
+%% 1c) Create the epsilon similarity graph
+[G] = epsilonSimGraph(epsilon,Pts, S);
+
+%% 1d) Create the adjacency matrix for the epsilon case
+W = S .* G;
+figure;
+gplotg(W,Pts(:,1:2))
+title('Adjacency matrix')
+%% 1e) Create the Laplacian matrix and implement spectral clustering
+[L,Diag] = CreateLapl(W);
+
+% \----------------------------/
+%     Your implementation
+% \----------------------------/
+
+% Cluster rows of eigenvector matrix of L corresponding to K smallest
+% eigennalues. Use kmeans for that.
+[D_spec,x_spec]     = kmeans_mod(Pts,K,n);
+
+%% 1f) Run K-means on input data
+% \----------------------------/
+%     Your implementation
+% \----------------------------/
+[D_kmeans,x_kmeans] = kmeans_mod(Pts,K,n);
+
+%% 1g) Visualize spectral and k-means clustering results
+figure;
+subplot(1,2,1)
+gplotmap(W,Pts,dummy_map)
+title('TODO: Plot the spectral clusters')
+subplot(1,2,2)
+gplotmap(W,Pts,dummy_map)
+title('TODO: Plot the K-means clusters')

mp4/Project_4_Maggioni_Claudio/src/CreateLapl.m

@@ -0,0 +1,18 @@
+function [L,Diag] = CreateLapl(W)
+% Create the Laplacian matrix of a graph
+% Input
+% W    : Adjacency matrix
+% Output
+% L    : Laplacian
+
+
+% Degree matrix
+Diag = zeros(size(W,1));
+for i = 1:size(W,1)
+   Diag(i,i) = sum(W(:,i)); 
+end
+
+% Construct the Laplacian
+L = Diag - W;
+
+end

mp4/Project_4_Maggioni_Claudio/src/epsilonSimGraph.m

@@ -0,0 +1,19 @@
+function [G] = epsilonSimGraph(epsilon, Pts, S)
+% Construct an epsilon similarity graph
+% Input
+% epsilon: size of neighborhood (calculate from Prim's Algorithm) 
+% Pts    : coordinate list of the sample 
+% 
+% Output
+% A      : the epsilon similarity matrix
+% 
+% USI, ICS, Lugano
+% Numerical Computing
+
+fprintf('----------------------------\n');
+fprintf('epsilon similarity graph\n');
+fprintf('----------------------------\n');
+
+G = S <= epsilon;
+
+end

mp4/Project_4_Maggioni_Claudio/src/gplotLength.m

@@ -0,0 +1,48 @@
+function [d] = gplotLength(a, x, linewidth)
+% function [d] = gplotLength(a, x,linewidth)
+%
+% plot with colors depending on lengths of
+% edges in x,
+% a: Adjacency
+% x: coords
+% linewidth: for ploting
+
+[ai,aj] = find(a);
+m = length(ai);
+d = double(a);
+
+for i = 1:m, 
+    d(ai(i),aj(i)) = norm(x(ai(i),:)-x(aj(i),:)); 
+end
+
+[di,dj,dv] = find(d);
+
+sdv = sort(dv);
+l = length(sdv);
+
+ran = [64:-1:0]/64;
+ran = ran.^2;
+ran = 1-ran;
+
+ll = round([1:l/65:l]);
+ranges = sdv(ll);
+
+clf; hold on;
+set(gca,'Color',[0 0 0])
+
+cm = jet;
+
+%cm = .3+.7*cm;
+
+for i = 1:64;
+ind = (dv <= ranges(i+1))&(dv>ranges(i));
+
+thisa = sparse(ai(ind),aj(ind),1);
+[gpx,gpy] = gplot(thisa,x);
+if (~isempty(gpx)),
+  ph(i) = plot(gpx,gpy,'-','Color',cm(i,:),'LineWidth',linewidth);
+  axis off
+end
+
+end
+

mp4/Project_4_Maggioni_Claudio/src/gplotg.m

@@ -0,0 +1,61 @@
+function handle = gplotg(A,xy,lc)
+% GPLOTG : Plot a "graph theoretic" graph.
+%
+%	handle = gplotg(A,xy) plots the graph specified by A and xy.
+%	A graph, G, is a set of nodes numbered from 1 to n,
+%	and a set of connections, or edges, between them.
+%	In order to plot G, two matrices are needed.
+%	The adjacency matrix, A, has a(i,j) nonzero if and
+%	only if node i is connected to node j.  The coordinates
+%	array, xy, is an n-by-2 or n-by-3  matrix with the position 
+%       for node i in the i-th row, xy(i,:) = [x(i) y(i)],
+%                               or  xy(i,:) = [x(i) y(i) z(i)].
+%	
+%	gplotg(A,xy,lc) uses line type and color instead of the 
+%	default, 'r-'.   For example, lc = 'g:'.  See PLOT.
+%
+%       Unlike gplot, gplotg sets 'erasemode' to 'none'
+%       so a graph won't redraw if something is drawn over it.
+%       Also, gplotg returns the handle to the plot.
+%       Also, gplotg can handle graphs with 3D coordinates.
+%	
+%	See also: SPY, TPLOT, SUBMESH, UNMESH.
+%
+%	John Gilbert, 1991.
+%	Modified 1-21-91, LS; 2-28-92, CBM; 6-14-01, JRG;
+%	Copyright (c) 1991-92 by the MathWorks, Inc.
+%       John Gilbert and Shanghua Teng, 1992-1993.
+% Copyright (c) 1990-1996 by Xerox Corporation.  All rights reserved.
+% HELP COPYRIGHT for complete copyright and licensing notice.
+
+if nargin < 3, 
+    lc = 'r-';
+end;
+
+[i,j] = find(A);
+[ignore, p] = sort(max(i,j));
+i = i(p);
+j = j(p);
+
+% Create a long, NaN-seperated list of line segments,
+% rather than individual segments.
+
+X = [ xy(i,1) xy(j,1) NaN*ones(size(i))]';
+Y = [ xy(i,2) xy(j,2) NaN*ones(size(i))]';
+X = X(:);
+Y = Y(:);
+
+if size(xy,2) == 2
+    h = plot (X, Y, lc, 'erasemode', 'none');
+else
+    set(gca,'drawmode','fast');
+    Z = [ xy(i,3) xy(j,3) NaN*ones(size(i))]';
+    Z = Z(:);
+    h = plot3 (X, Y, Z, lc, 'erasemode', 'none');
+end;
+
+axis equal;
+axis off;
+if nargout >= 1
+    handle = h;
+end;

mp4/Project_4_Maggioni_Claudio/src/gplotmap.m

@@ -0,0 +1,86 @@
+function handle = gplotmap(A,xy,map,pcolor,ecolor)
+% GPLOTMAP : Plot a partitioned graph in 2 or 3 dimensions.
+%
+%	gplotmap(A,xy,map) plots the n-vertex graph specified 
+%       by the n by n adjacency (or Laplacian) matrix A 
+%       and the n by 2 or 3 matrix of coordinates xy.
+%       Argument map is an n-vector of part numbers, which
+%       assigns each vertex to a part.
+%
+%       By default, edges that join different parts are omitted, and
+%       the picture shows each part in a different color.  The call
+%
+%       gplotmap(A,xy,map,pcolor,ecolor)
+%
+%       uses color "pcolor" for the parts and "ecolor" for the
+%       edges between the parts.  If "pcolor" has multiple rows,
+%       each part gets the color of one row (in cyclic order).
+%
+% See also GPLOTPART.
+%
+% John Gilbert  2 August 1994
+% Copyright (c) 1990-1996 by Xerox Corporation.  All rights reserved.
+% HELP COPYRIGHT for complete copyright and licensing notice.
+
+[n,n] = size(A);
+if nargin < 3, 
+    map = zeros(1,n);
+end;
+if nargin < 4 
+    pcolor = [];
+end;
+if nargin < 5,
+    ecolor = [];
+end;
+
+parts = setfilter(map);
+nparts = length(parts);
+if length(pcolor) == 0,
+    pcolor = hsv(nparts);
+    pcolor = pcolor(randperm(nparts),:);
+end;
+
+% clf reset
+% colordef(gcf,'black')
+if size(xy,2) == 2
+    axis([min(xy(:,1)) max(xy(:,1)) min(xy(:,2)) max(xy(:,2))]);
+else
+    axis([min(xy(:,1)) max(xy(:,1)) min(xy(:,2)) max(xy(:,2)) ...
+        min(xy(:,3)) max(xy(:,3))]); 
+end;
+
+axis equal;
+axis off;
+lighting gouraud;
+hold on
+
+% Count and plot the separating edges.
+[i,j] = find(A);
+f = find(map(i) > map(j));
+if length(f) 
+%     xlabel( [int2str(length(f)) ' cut edges on ' int2str(nparts) ' partitions'],'visible','on');
+    if length(ecolor)
+        cut = sparse(i(f),j(f),1,n,n);
+        set(gplotg(cut,xy,'-'),'color',ecolor);
+    end;
+else
+%     xlabel('0 cut edges','visible','on');
+end;
+
+% Plot each piece.
+ncolor = 1;
+for partnumber = parts;
+    c = pcolor(ncolor,:);
+    ncolor = rem(ncolor, size(pcolor,1)) + 1;
+    part = find(map == partnumber);
+    set(gplotg(A(part,part),xy(part,:),'-'),'color',c);
+    if n < 500,
+        if size(xy,2) == 2
+            set(plot(xy(part,1),xy(part,2),'o'),'color',c);
+        else
+            set(plot3(xy(part,1),xy(part,2),xy(part,3),'o'),'color',c);
+        end;
+    end;
+end;
+
+hold off

mp4/Project_4_Maggioni_Claudio/src/isConnected.m

@@ -0,0 +1,70 @@
+%##################################################################
+% Determine if a graph is connected
+% Idea by Ed Scheinerman, circa 2006, source: http://www.ams.jhu.edu/~ers/matgraph/
+%                                     routine: matgraph/@graph/isconnected.m
+%
+% INPUTS: adjacency matrix, nxn
+% OUTPUTS: Boolean variable, 0 or 1
+%
+% Note: This function works only for undirected graphs.
+% GB: last updated, Sep 23 2012
+%##################################################################
+
+function S = isConnected(adj)
+
+if not(isempty(find(sum(adj)==0))); S = false; return; end
+
+n = length(adj);
+x = [1; zeros(n-1,1)]; % [1,0,...0] nx1 vector 
+
+while 1
+     y = x;
+     x = adj*x + x;
+     x = x>0;
+     
+     if x==y; break; end
+
+end
+
+S = true;
+if sum(x)<n; S = false; 
+
+end
+
+end
+
+
+% Alternative 1 ==========================================================
+% If the algebraic connectivity is > 0 then the graph is connected
+% a=algebraic_connectivity(adj);
+% S = false; if a>0; S = true; end
+
+% Alternative 2 ==========================================================
+% Uses the fact that multiplying the adj matrix to itself k times give the
+% number of ways to get from i to j in k steps. If the end of the
+% multiplication in the sum of all matrices there are 0 entries then the
+% graph is disconnected. Computationally intensive, but can be sped up by
+% the fact that in practice the diameter is very short compared to n, so it
+% will terminate in order of log(n)? steps.
+% function S=isconnected(el):
+%     
+%     S=false;
+%     
+%     adj=edgeL2adj(el);
+%     n=numnodes(adj); % number of nodes
+%     adjn=zeros(n);
+% 
+%     adji=adj;
+%     for i=1:n
+%         adjn=adjn+adji;
+%         adji=adji*adj;
+% 
+%         if length(find(adjn==0))==0
+%             S=true;
+%             return
+%         end
+%     end
+
+% Alternative 3 ==========================================================
+% Find all connected components, if their number is 1, the graph is
+% connected. Use find_conn_comp(adj).

mp4/Project_4_Maggioni_Claudio/src/kNNSimGraph.m

@@ -0,0 +1,28 @@
+function [G] = kNNSimGraph(Pts)    
+% Construct a k-nearest neighbors similarity graph
+% Input
+% k      : # of neighbors
+% Pts    : coordinate list of the sample 
+% 
+% Output
+% G      : the kNN similarity matrix
+
+
+n = length(Pts(:,1));
+kn = ceil(2*log(n));
+
+
+fprintf('kNN similarity graph\n');
+    for i = 1:n
+        s = repmat(Pts(i,:),n,1);
+        d = Pts - s;
+%         e = diag(d*d');
+        e = sum(d.^2,2);
+        [val,ind] = sort(e);
+        ind(1) = [];
+        nbrs = ind(1:kn);
+        G(i,nbrs) = 1;
+        G(nbrs,i) = 1;
+    end
+    
+end

mp4/Project_4_Maggioni_Claudio/src/kmeans_mod.m

@@ -0,0 +1,57 @@
+function [D,x] = kmeans_mod(Y,K,n)
+% Function to perform k-means clustering
+% --------------------------------------
+% Input
+% Y: data matrix
+% K: # of clusters
+% n: number of points
+% Output
+% D: centroids
+% x: assignment(indicator) vector
+
+
+cand = unique(round(10000*Y)/10000,'rows');
+c = datasample(cand,K,'Replace',false); % pick K random input pts as initial centroids
+D = c;
+x = zeros(size(Y,1),1);
+D_old = inf*ones(size(D));
+count = 1;
+while norm(D - D_old) > 0.00001 & count < 500
+    D_old = D;
+    % Assign points to clusters
+    for i = 1:n
+        min = inf;
+        for k = 1:K
+            dist = norm(Y(i,:) - D(k,:));
+            if dist < min
+                min = dist;
+                x(i) = k;
+            end
+        end
+    end
+    
+    % Update centroid locations
+    for k = 1:K
+        ind = find(x == k);
+        if length(ind) == 1
+            D(k,:) = Y(ind,:);
+        else
+            D(k,:) = mean(Y(ind,:));
+        end
+    end
+    
+    count = count + 1;
+end
+
+for i = 1:n
+    min = inf;
+    for k = 1:K
+        dist = norm(Y(i,:) - D(k,:));
+        if dist < min
+            min = dist;
+            x(i) = k;
+        end
+    end
+end
+end
+

mp4/Project_4_Maggioni_Claudio/src/minSpanTree.m

@@ -0,0 +1,51 @@
+function A_sp = minSpanTree(A)
+%MINSPANTREE Prim's minimal spanning tree algorithm
+
+% @input A, NxN adjacency matrix
+% @output A_sp, NxN adjacency matrix of minimum spanning tree
+
+% Prim's alg idea:
+%  start at any node, find closest neighbor and mark edges
+%  for all remaining nodes, find closest to previous cluster, mark edge
+%  continue until no nodes remain
+%
+% INPUTS: graph defined by adjacency matrix, nxn
+% OUTPUTS: matrix specifying minimum spanning tree (subgraph), nxn
+%
+% Other routines used: isConnected.m
+% Updated: cleaned fprintf
+
+% IB Last updated: 3/24/14
+
+
+A_sp = [];
+
+% check if graph is connected:
+if not(isConnected(A)); fprintf('This graph is not connected. No spanning tree exists.\n'); return; end
+
+n = length(A); % number of nodes
+A_sp = sparse(n,n);   % initialize tree
+
+A(A==0)=inf; % set all zeros in the matrix to inf
+
+conn_nodes = 1;        % nodes part of the min-span-tree
+rem_nodes = [2:n];     % remaining nodes
+
+while ~isempty(rem_nodes)
+  
+  [minlink]=min(min(A(conn_nodes,rem_nodes)));
+  ind=find(A(conn_nodes,rem_nodes)==minlink);
+
+  [ind_i,ind_j] = ind2sub([length(conn_nodes),length(rem_nodes)],ind(1));
+
+  i=conn_nodes(ind_i); j=rem_nodes(ind_j); % gets back to adj indices
+  A_sp(i,j)=A(i,j); A_sp(j,i)=A(j,i);
+  conn_nodes = [conn_nodes j]; %#ok<AGROW>
+  rem_nodes = setdiff(rem_nodes,j);
+  
+end
+
+
+
+
+end

mp4/Project_4_Maggioni_Claudio/src/setfilter.m

@@ -0,0 +1,20 @@
+function B = setfilter(A)
+% SETFILTER : Sort and remove duplicates.
+%
+% B = setfilter(A) returns a row vector with one occurrence
+%     of each different element of A, in sorted order.
+%
+% Based on code by Asko Huuskonen, Finnish Meteorological Institute,
+% Asko.Huuskonen@fmi.fi, posted to comp.soft-sys.matlab.
+% This version by John Gilbert, Xerox PARC, 27 Sep 1993
+% Copyright (c) 1990-1996 by Xerox Corporation.  All rights reserved.
+% HELP COPYRIGHT for complete copyright and licensing notice.
+
+if length(A) == 0
+    B = []; 
+    return;
+end;
+
+B = sort(A(:));
+B(find(diff(B)==0)) = [];
+B = B.';

mp4/Project_4_Maggioni_Claudio/src/similarityfunc.m

@@ -0,0 +1,19 @@
+function [S] = similarityfunc(Pts, sigma)
+% Create the similarity matrix S from the coordinate list of the input
+% points
+% dimosthenis.pasadakis@usi.ch
+% ICS, USI.
+
+if nargin < 2
+    % Choose \sigma ~ 2*log(n)
+    n = length(Pts(:,1));
+    sigma = log(n);
+end
+
+fprintf('----------------------------\n');
+fprintf('Gaussian similarity function\n');
+fprintf('----------------------------\n');
+S = squareform(pdist(Pts));
+S = exp(-S.^2 ./ (2*sigma^2));
+
+end

mp4/assignment.sty

@@ -0,0 +1,95 @@
+\usepackage{ifthen}
+\usepackage[utf8]{inputenc}
+\usepackage{graphics}
+\usepackage{graphicx}
+\usepackage{hyperref}
+
+\pagestyle{plain}
+\voffset -5mm
+\oddsidemargin  0mm
+\evensidemargin -11mm
+\marginparwidth 2cm
+\marginparsep 0pt
+\topmargin 0mm
+\headheight 0pt
+\headsep 0pt
+\topskip 0pt        
+\textheight 255mm
+\textwidth 165mm
+
+\newcommand{\duedate} {}
+\newcommand{\setduedate}[1]{%
+\renewcommand\duedate {Due date:~ #1}}
+\newcommand\isassignment {false}
+\newcommand{\setassignment}{\renewcommand\isassignment {true}}
+\newcommand{\ifassignment}[1]{\ifthenelse{\boolean{\isassignment}}{#1}{}}
+\newcommand{\ifnotassignment}[1]{\ifthenelse{\boolean{\isassignment}}{}{#1}}
+
+\newcommand{\assignmentpolicy}{
+\begin{table}[h]
+\begin{center}
+\scalebox{0.8} {%
+\begin{tabular}{|p{0.02cm}p{16cm}|}
+\hline
+&\\
+\multicolumn{2}{|c|}{\Large\textbf{Numerical Computing 2020 ---  Submission Instructions}}\\
+\multicolumn{2}{|c|}{\large\textbf{(Please, notice that following instructions are mandatory: }}\\
+\multicolumn{2}{|c|}{\large\textbf{submissions that don't comply with, won't be considered)}}\\
+&\\
+\textbullet & Assignments must be submitted to \href{https://www.icorsi.ch/course/view.php?id=10018}{iCorsi} (i.e. in electronic format).\\
+\textbullet & Provide both executable package and sources (e.g. C/C++ files, Matlab). 
+If you are using libraries, please add them in the file. Sources must be organized in directories called:\\
+\multicolumn{2}{|c|}{\textit{Project\_number\_lastname\_firstname}}\\
+& and  the  file must be called:\\
+\multicolumn{2}{|c|}{\textit{project\_number\_lastname\_firstname.zip}}\\
+\multicolumn{2}{|c|}{\textit{project\_number\_lastname\_firstname.pdf}}\\
+\textbullet &  The TAs will grade your project by reviewing your project write-up, and looking at the implementation 
+                 you attempted, and benchmarking your code's performance.\\
+
+\textbullet & You are allowed to discuss all questions with anyone you like; however: (i) your submission must list anyone you discussed problems with and (ii) you must write up your submission independently.\\
+\hline
+\end{tabular}
+}
+\end{center}
+\end{table}
+}
+\newcommand{\punkte}[1]{\hspace{1ex}\emph{\mdseries\hfill(#1~\ifcase#1{Points}\or{Points}\else{Points}\fi)}}
+
+
+\newcommand\serieheader[6]{
+\thispagestyle{empty}%
+\begin{flushleft}
+\includegraphics[width=0.4\textwidth]{usi_inf}
+\end{flushleft}
+  \noindent%
+  {\large\ignorespaces{\textbf{#1}}\hspace{\fill}\ignorespaces{ \textbf{#2}}}\\ \\%
+  {\large\ignorespaces #3 \hspace{\fill}\ignorespaces #4}\\
+  \noindent%
+  \bigskip
+  \hrule\par\bigskip\noindent%
+  \bigskip {\ignorespaces {\Large{\textbf{#5}}}
+  \hspace{\fill}\ignorespaces \large \ifthenelse{\boolean{\isassignment}}{\duedate}{#6}}
+  \hrule\par\bigskip\noindent%  \linebreak
+ }
+
+\makeatletter
+\def\enumerateMod{\ifnum \@enumdepth >3 \@toodeep\else
+      \advance\@enumdepth \@ne
+      \edef\@enumctr{enum\romannumeral\the\@enumdepth}\list
+      {\csname label\@enumctr\endcsname}{\usecounter
+        {\@enumctr}%%%? the following differs from "enumerate"
+	\topsep0pt%
+	\partopsep0pt%
+	\itemsep0pt%
+	\def\makelabel##1{\hss\llap{##1}}}\fi}
+\let\endenumerateMod =\endlist
+\makeatother
+
+
+
+
+\usepackage{textcomp}
+
+
+
+