clear all;
close all;
clc;

%Caricamento delle matrici O/D e della mappatura
load 'IMQ2013_zonecampionarie';
fprintf('File MAT di riferimento:  IMQ2013_zonecampionarie.mat\n');
 
%Ottenimento delle informazioni sul criterio di aggregazione
[empty,file] = xlsread('AGGREGAZIONE.csv','A2:C210');
zoneORD = string(file(:,1)); zoneORD = [zoneORD;'Esterno'];
PROV = string(file(:,2)); 
MACR = string(file(:,3)); 
CRIT = input('\nInserire criterio di aggregazione (1-per province, 2-per macro-aree): ');

%Aggregazione delle matrici
if CRIT==1;
    NEWzone = unique(PROV);
    PROV = [PROV;'EXT']; NEWzone = [NEWzone;'EXT'];
    [AGGR_ResPIEmPrIGt] = aggregazione(ResPIEmPrIGt,NEWzone,PROV,zoneORD,zoneA);
    [AGGR_ResPIEmPrIGsr] = aggregazione(ResPIEmPrIGsr,NEWzone,PROV,zoneORD,zoneA);
    [AGGR_ResTOmPrIGt] = aggregazione(ResTOmPrIGt,NEWzone,PROV,zoneORD,zoneA);
    [AGGR_ResTOmPrIGsr] = aggregazione(ResTOmPrIGsr,NEWzone,PROV,zoneORD,zoneA);
    [AGGR_ResPIEmPrOPt] = aggregazione(ResPIEmPrOPt,NEWzone,PROV,zoneORD,zoneA);
    [AGGR_ResTOmPrOPt] = aggregazione(ResTOmPrOPt,NEWzone,PROV,zoneORD,zoneA);
    [AGGR_ResPIEmPuIGt] = aggregazione(ResPIEmPuIGt,NEWzone,PROV,zoneORD,zoneA);
    [AGGR_ResPIEmPuIGsr] = aggregazione(ResPIEmPuIGsr,NEWzone,PROV,zoneORD,zoneA);
    [AGGR_ResTOmPuIGt] = aggregazione(ResTOmPuIGt,NEWzone,PROV,zoneORD,zoneA);
    [AGGR_ResTOmPuIGsr] = aggregazione(ResTOmPuIGsr,NEWzone,PROV,zoneORD,zoneA);
    [AGGR_ResPIEmPuOPt] = aggregazione(ResPIEmPuOPt,NEWzone,PROV,zoneORD,zoneA);
    [AGGR_ResTOmPuOPt] = aggregazione(ResTOmPuOPt,NEWzone,PROV,zoneORD,zoneA);
else if CRIT==2;
        NEWzone = unique(MACR);
        MACR = [MACR;'EXT']; NEWzone = [NEWzone;'EXT'];
        [AGGR_ResPIEmPrIGt] = aggregazione(ResPIEmPrIGt,NEWzone,MACR,zoneORD,zoneA);
        [AGGR_ResPIEmPrIGsr] = aggregazione(ResPIEmPrIGsr,NEWzone,MACR,zoneORD,zoneA);
        [AGGR_ResTOmPrIGt] = aggregazione(ResTOmPrIGt,NEWzone,MACR,zoneORD,zoneA);
        [AGGR_ResTOmPrIGsr] = aggregazione(ResTOmPrIGsr,NEWzone,MACR,zoneORD,zoneA);
        [AGGR_ResPIEmPrOPt] = aggregazione(ResPIEmPrOPt,NEWzone,MACR,zoneORD,zoneA);
        [AGGR_ResTOmPrOPt] = aggregazione(ResTOmPrOPt,NEWzone,MACR,zoneORD,zoneA);
        [AGGR_ResPIEmPuIGt] = aggregazione(ResPIEmPuIGt,NEWzone,MACR,zoneORD,zoneA);
        [AGGR_ResPIEmPuIGsr] = aggregazione(ResPIEmPuIGsr,NEWzone,MACR,zoneORD,zoneA);
        [AGGR_ResTOmPuIGt] = aggregazione(ResTOmPuIGt,NEWzone,MACR,zoneORD,zoneA);
        [AGGR_ResTOmPuIGsr] = aggregazione(ResTOmPuIGsr,NEWzone,MACR,zoneORD,zoneA);
        [AGGR_ResPIEmPuOPt] = aggregazione(ResPIEmPuOPt,NEWzone,MACR,zoneORD,zoneA);
        [AGGR_ResTOmPuOPt] = aggregazione(ResTOmPuOPt,NEWzone,MACR,zoneORD,zoneA);
    end;
end;

%Salvataggio delle nuove matrici
save IMQ2013_aggregazione NEWzone AGGR_ResPIEmPrIGt AGGR_ResPIEmPrIGsr AGGR_ResTOmPrIGt AGGR_ResTOmPrIGsr AGGR_ResPIEmPrOPt AGGR_ResTOmPrOPt AGGR_ResPIEmPuIGt AGGR_ResPIEmPuIGsr AGGR_ResTOmPuIGt AGGR_ResTOmPuIGsr AGGR_ResPIEmPuOPt AGGR_ResTOmPuOPt;
fprintf('\nFILE MAT CREATO:  IMQ2013_aggregazione.mat\n');

%Definizione della funzione di aggregazione
function [newMATR] = aggregazione(matrORIG,Nzone,criterio,zORDIN,zORIG);
newMATR = zeros(length(Nzone),length(Nzone));
for i=1:209;
    indexO = find(zORDIN==zORIG(i));
    iO = find(Nzone==criterio(indexO));
    for j=1:209;
        indexD = find(zORDIN==zORIG(j));
        iD = find(Nzone==criterio(indexD));
        newMATR(iO,iD) = newMATR(iO,iD)+matrORIG(i,j);
    end;
end;
end