Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 0a8babdc authored by hhakim's avatar hhakim
Browse files

Add FGFT error benchmark scripts (python and matlab) and a graph generation...

Add FGFT error benchmark scripts (python and matlab) and a graph generation script using GSP matlab toolbox.

The goal of these scripts is to validate the C++ port of the FGFT algorithms and their python wrapper against the reference Matlab prototypes from FAuST 1.03.
The two FGFT algorithms are tested : PALM FGFT and Givens FGFT.
parent 1e330111
No related branches found
No related tags found
No related merge requests found
import matfaust.Faust
fpath = mfilename('fullpath');
[path, ~, ~ ] = fileparts(fpath);
fs = filesep;
data_path = [ path fs '..' fs '..' fs '..' fs 'misc' fs 'data' fs 'mat' fs 'Laps_U_Ds-6_graph_types-dims_128_to_1024.mat'];
if(~ exist(data_path))
error([data_path ' not found. Please run the data generation script misc/test/src/Matlab/get_graph_test_UDLap_matrices.m'])
end
load(data_path)
num_graphs = length(Laps);
hier_errs = {};
givens_errs = {};
par_givens_errs = {};
U_hier_errs = {};
U_givens_errs = {};
U_par_givens_errs = {};
num_graphs = 60; % only first 60 graph Laps which are of 6 types (random sensor, etc.)
for i=1:num_graphs
disp(['benchmarking graph. Laplacian ' num2str(i) '/' num2str(num_graphs)])
U = Us{i};
D = Ds{i};
Lap = Laps{i};
dim = size(Lap, 1);
nfacts = round(log2(dim))-3; % Desired number of factors
over_sp = 1.5; % Sparsity overhead
dec_fact = 0.5; % Decrease of the residual sparsity
params.data = U;
params.Lap = Lap;
params.init_D = D;
params.nfacts = nfacts;
params.fact_side = 0;
params.cons = cell(2,nfacts-1);
for j=1:nfacts-1
params.cons{1,j} = {'sp',dec_fact^j*dim^2*over_sp,dim,dim};
params.cons{2,j} = {'sp',2*dim*over_sp,dim,dim};
end
% Number of iterations
params.niter1 = 50;
params.niter2 = 100;
params.verbose=0;
[lambda, facts, errors] = hierarchical_fact(params);
complexity_global=0;
for j=1:nfacts
complexity_global = complexity_global + nnz(facts{j});
end
F = Faust(facts);
RC_PALM = complexity_global/nnz(U);
[lambda, facts, Dhat, errors] = hierarchical_fact_FFT(params);
facts{1} = lambda*facts{1};
Uhat_PALM = Faust(facts);
hier_err = norm(Uhat_PALM*Dhat*Uhat_PALM' - Lap, 'fro')/norm(Lap, 'fro');
hier_errs = [ hier_errs hier_err];
U_hier_errs = [ U_hier_errs norm(Uhat_PALM - U, 'fro')/norm(U, 'fro')];
J=round(RC_PALM*nnz(U)/4);
[facts_givens,Dhat,err,L,choices] = diagonalization_givens(Lap,J);
Uhat_givens = Faust(facts_givens);
givens_err = norm(Uhat_givens*full(Dhat)*Uhat_givens' - Lap, 'fro')/norm(Lap, 'fro');
givens_errs = [ givens_errs givens_err];
U_givens_errs = [ U_givens_errs norm(Uhat_givens - U, 'fro')/norm(U, 'fro')];
[facts_givens_parall,Dhat,err,L,coord_choices] = diagonalization_givens_parall(Lap,J,dim/2);
Uhat_givens_par = Faust(facts_givens_parall);
par_givens_err = norm(Uhat_givens_par*full(Dhat)*Uhat_givens_par' - Lap, 'fro')/norm(Lap, 'fro');
par_givens_errs = [ par_givens_errs par_givens_err];
U_par_givens_errs = [ U_par_givens_errs norm(Uhat_givens_par - U, 'fro')/norm(U, 'fro')];
end
save('benchmark_Lap_diag_output.mat', 'hier_errs', 'givens_errs', 'par_givens_errs', 'U_hier_errs', 'U_givens_errs', 'U_par_givens_errs')
hier_errs
givens_errs
par_givens_errs
U_hier_errs
U_givens_errs
U_par_givens_errs
% code inspired from FAuST 1.03 /experiment_comparison_PALM_givens.m to extract graph tested matrices
% necessitates the GSP toolbox to generates graphs/Laplacians
DIMtab = [128,256,512,1024]; % Dimension of the graphs (number of nodes)
Nrun=10;
resTab = zeros(4, numel(DIMtab));
graph_types = {'erdos_renyi', 'community', 'sensor', 'path', 'random_ring', 'ring'}
Laps = {}
Us = {}
Ds = {}
for k = 1:numel(DIMtab)
DIM = DIMtab(k);
%disp(DIM)
for n=1:Nrun
%disp(n)
for l = 1:6 % l == 3
disp(['k=' num2str(k) ])%, n=' num2str(n)]) %)', l=' num2str(l)])
%%%
%% ==================== Graph generation =========================
%%%
if l==1
%disp(['l = ' num2str(l)])
G = gsp_erdos_renyi(DIM,0.1);
elseif l==2
%disp(['l = ' num2str(l)])
G = gsp_community(DIM);
elseif l==3
%disp(['l = ' num2str(l)])
G = gsp_sensor(DIM);
elseif l==4
%disp(['l = ' num2str(l)])
G = gsp_path(DIM);
elseif l==5
%disp(['l = ' num2str(l)])
G = gsp_random_ring(DIM);
elseif l==6
%disp(['l = ' num2str(l)])
G = gsp_ring(DIM);
end
Lap = full(G.L);
%save(['Laplacian_' num2str(DIMtab(k)) '_' graph_types{l}], 'Lap')
Laps = [ Laps Lap ]
% G = gsp_compute_fourier_basis(G);
%% U = full(G.U);
%% D = diag(G.e);
[U,D]=eig(Lap);
Us = [ Us U ]
Ds = [ Ds D ]
end
end
end
save('Laps_U_Ds-6_graph_types-dims_128_to_1024.mat', 'Laps', 'Us', 'Ds')
from scipy.io import loadmat
from sys import argv
from os.path import dirname, sep, join, exists
from numpy import diag, copy, log2, count_nonzero, size, loadtxt, savetxt
import numpy, re
from numpy.linalg import norm
from pyfaust.factparams import *
from pyfaust import *
from pylab import *
if __name__ == '__main__':
data_file = join(dirname(argv[0]),
'../../../data/mat/Laps_U_Ds-6_graph_types-dims_128_to_1024.mat')
if(not exists(data_file)):
raise Exception(data_file+" not found. You need to run"
" misc/test/src/Matlab/get_graph_test_UDLap_matrices.m"
" for generation.")
matfile = loadmat(data_file)
Laps = matfile['Laps'][0]
Ds = matfile['Ds'][0]
Us = matfile['Us'][0]
start_i=0
num_laps = Laps.shape[0]
num_laps = 45
hier_fgft_errs = empty(num_laps)
givens_errs = empty(num_laps)
par_givens_errs = empty(num_laps)
U_hier_fgft_errs = empty(num_laps)
U_givens_errs = empty(num_laps)
U_par_givens_errs = empty(num_laps)
if exists('benchmark_pyfaust_fgft.txt'):
saved_lap_errs = loadtxt('benchmark_pyfaust_fgft.txt')
len_saved = saved_lap_errs[0].shape[0]
givens_errs[:len_saved], par_givens_errs[:len_saved],\
hier_fgft_errs[:len_saved] = saved_lap_errs
saved_U_errs = loadtxt('benchmark_pyfaust_fgft_U.txt')
U_givens_errs[:len_saved], U_par_givens_errs[:len_saved], \
U_hier_fgft_errs[:len_saved] = \
saved_U_errs
start_i = len_saved
if(saved_U_errs[0].shape[0] != saved_lap_errs[0].shape[0]):
raise Exception("Error files not consistent, delete and recompute")
for i in range(start_i, num_laps):
print('Lap ',i,'/',num_laps)
U,D,Lap = Us[i],Ds[i],Laps[i].astype(numpy.float)
dim = Lap.shape[0]
nfacts = int(round(log2(dim))-3)
over_sp = 1.5 # sparsity overhead
dec_fact = .5 # decrease of the residum sparsity
fact_cons, res_cons = [], []
for j in range(1, nfacts):
fact_cons += [ ConstraintInt('sp', dim, dim,
min(int(round(dec_fact**j*dim**2*over_sp)),size(Lap)))
]
res_cons += [ ConstraintInt('sp', dim, dim,
min(int(round(2*dim*over_sp)),size(Lap))) ]
params = ParamsHierarchicalFact(fact_cons, res_cons,
StoppingCriterion(num_its=50),
StoppingCriterion(num_its=100),
step_size=1.0000e-06,
constant_step_size=False,
init_lambda=1.0,
is_fact_side_left=False)
F = FaustFactory.fact_hierarchical(U,params)
print(F)
complexity_global = F.nnz_sum()
rc_palm = complexity_global / count_nonzero(U)
diag_init_D = diag(D)
diag_init_D = numpy.copy(diag_init_D)
F, Dhat = FaustFactory.fgft_palm(U, Lap,
params,
diag_init_D)
hier_fgft_err = norm((F.todense()*diag(Dhat))*F.T.todense()-Lap,"fro")/norm(Lap,"fro")
hier_fgft_errs[i] = hier_fgft_err
U_hier_fgft_errs[i] = (F-U).norm("fro")/norm(U,"fro")
#nfacts = complexity_global/(2*dim)
#J = round(nfacts*(dim/2))
J = round(complexity_global/4)
F, Dhat = FaustFactory.fgft_givens(Lap, J, 0)
givens_err = norm((F*Dhat.todense())*F.T.todense()-Lap,'fro')/norm(Lap,'fro')
givens_errs[i] = givens_err
U_givens_errs[i] = (F-U).norm("fro")/norm(U,"fro")
F, Dhat = FaustFactory.fgft_givens(Lap, J, int(dim/2))
print("J=", J)
print("nnz_sum FGFT givens parallel=", F.nnz_sum(), "num of facts:",
F.get_num_factors())
par_givens_err = norm((F*Dhat.todense())*F.T.todense()-Lap,'fro')/norm(Lap,'fro')
par_givens_errs[i] = par_givens_err
U_par_givens_errs[i] = (F-U).norm("fro")/norm(U,"fro")
print("err:", norm(Lap * F.toarray() - F.toarray()*Dhat)/norm(Lap))
savetxt('benchmark_pyfaust_fgft.txt', np.array([givens_errs[:i+1],
par_givens_errs[:i+1],
hier_fgft_errs[:i+1]]))
savetxt('benchmark_pyfaust_fgft_U.txt', np.array([U_givens_errs[:i+1],
U_par_givens_errs[:i+1],
U_hier_fgft_errs[:i+1]]))
matlab_matfile = None
if(len(argv) > 1 and exists(argv[1]) and re.match(".*.mat", argv[1])):
matlab_matfile = loadmat(argv[1])
matlab_hier_fgft_errs = [matlab_matfile['hier_errs'][0,i][0,0] for i in
range(matlab_matfile['hier_errs'].shape[1]) ]
matlab_givens_errs = [matlab_matfile['givens_errs'][0,i][0,0] for i in
range(matlab_matfile['givens_errs'].shape[1]) ]
matlab_par_givens_errs = [matlab_matfile['par_givens_errs'][0,i][0,0] for i in
range(matlab_matfile['par_givens_errs'].shape[1]) ]
print('hier_fgft_errs = ', hier_fgft_errs)
print('givens_errs = ', givens_errs)
print('par_givens_errs = ', par_givens_errs)
print('U_hier_fgft_errs = ', U_hier_fgft_errs)
print('U_givens_errs = ', U_givens_errs)
print('U_par_givens_errs = ', U_par_givens_errs)
plt.rcParams['figure.figsize'] = [18.0, 12]
lap_indices = arange(num_laps)
grid(True, which='both', axis='both')
xticks(lap_indices)
scatter(lap_indices, givens_errs, c='b', marker='o', label='pyfaust_givens_errs')
scatter(lap_indices, par_givens_errs, c='r', marker='o', label='pyfaust_par_givens_errs')
scatter(lap_indices, hier_fgft_errs, c='y', marker='o', label='pyfaust_hier_fgft_errs')
if matlab_matfile:
scatter(lap_indices, matlab_hier_fgft_errs[:len(lap_indices)],
label='matlab_hier_fgft_errs', c='y', marker='+', s=315)
scatter(lap_indices, matlab_givens_errs[:len(lap_indices)],
label='matlab_givens_errs', c='b', marker='+', s=315)
scatter(lap_indices, matlab_par_givens_errs[:len(lap_indices)],
label='matlab_par_givens_errs', c='r', marker='+', s=315)
ylabel("Relative error on Laplacian (reconstructed with FGFT)")
xlabel("Laplacians\n i-th Lap. is for"
" erdos-renyi if i==0 (mod 6)\ncommunity if i==1 (mod 6)\nsensor if i=="
"2 (mod 6)\npath if i==3 (mod 6)\nrandom ring if i==4 (mod6), ring if i=="
"5 (mod6)")
title("Error benchmark on "+str(num_laps)+" Laplacians (128x128)")
legend()
savefig('benchmark_Lap_diag_pyfaust_Laplacian_figure.png')
figure(2)
grid(True, which='both', axis='both')
if matlab_matfile:
matlab_U_hier_fgft_errs = [matlab_matfile['U_hier_errs'][0,i][0,0] for i in
range(matlab_matfile['U_hier_errs'].shape[1]) ]
matlab_U_givens_errs = [matlab_matfile['U_givens_errs'][0,i][0,0] for i in
range(matlab_matfile['U_givens_errs'].shape[1]) ]
matlab_U_par_givens_errs = [matlab_matfile['U_par_givens_errs'][0,i][0,0] for i in
range(matlab_matfile['U_par_givens_errs'].shape[1]) ]
scatter(lap_indices, matlab_U_hier_fgft_errs[:len(lap_indices)],
label='matlab_hier_fgft_errs', c='y', marker='+', s=315)
scatter(lap_indices, matlab_U_givens_errs[:len(lap_indices)],
label='matlab_givens_errs', c='b', marker='+', s=315)
scatter(lap_indices, matlab_U_par_givens_errs[:len(lap_indices)],
label='matlab_par_givens_errs', c='r', marker='+', s=315)
xticks(lap_indices)
scatter(lap_indices, U_givens_errs, c='b', marker='o', label='pyfaust_givens_errs')
scatter(lap_indices, U_par_givens_errs, c='r', marker='o', label='pyfaust_par_givens_errs')
scatter(lap_indices, U_hier_fgft_errs, c='y', marker='o', label='pyfaust_hier_fgft_errs')
ylabel("Relative error of FGFT versus U (Fourier matrix)")
xlabel("Fourier matrices")
title("Error benchmark on "+str(num_laps)+" Graph Fourier matrices (128x128)")
legend()
if matlab_matfile:
print('matlab_hier_fgft_errs = ', matlab_hier_fgft_errs)
print('matlab_givens_errs = ', matlab_givens_errs)
print('matlab_par_givens_errs = ', matlab_par_givens_errs)
print('matlab_U_hier_fgft_errs = ', matlab_U_hier_fgft_errs)
print('matlab_U_givens_errs = ', matlab_U_givens_errs)
print('matlab_U_par_givens_errs = ', matlab_U_par_givens_errs)
savefig('benchmark_Lap_diag_pyfaust_Fourier_figure.png')
show()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment