% IANTA15: Programa MATLAB para "Proyecto IANTA"
%
% Representación temporal y espectral del sonido de anuros
%
% 15: Clasificación con número variable de frames
%       Se toman los parámetros principales de varios frames
%

function IANTA15

clear all;
clear global all;
clc;

global PATRON;
global carpeta; %Carpeta de trabajo de los sonidos
global ParUtil; %Parámetros utilizados
global NumPar NumFrames;

%-------------- Variables generales
[FicheroPatron,roiPatron,FicheroResumenPatrones]=DefinePatron;
DefineVariablesGenerales;
carpeta='sonidos\Procesados\';
%carpeta='sonidos\kk\';

OrdPar=[6,1,8,10,7,12,5,3,18,4,17,16,14,13,11,15,9,2];
NumPar=5; %Número de parámetros utilizados
NumFrames=5; %Número de frames utilizados
ParUtil=OrdPar(1:NumPar);
FicheroResumen=[carpeta,'ResumenSW-',num2str(NumPar)];

%{
[matparx,ClasePatronNumx,NomPar]=GeneraDatosPatrones;
nobs=size(matparx,1);
matpar=zeros(nobs-NumFrames+1,NumFrames*NumPar);
for k=NumFrames:nobs
    matx=[];
    k0=k-NumFrames+1;
    for k1=k0:k
        x=matparx(k1,ParUtil);
        matx=[matx,x];
    end
    matpar(k0,:)=matx;
end
ClasePatronNum=ClasePatronNumx(NumFrames:end);
%}

%{
%-------------- Calcula la clase de sonido real de cada fichero
[NombreFichero,ClaseSonidoFichero]=ClasificaCarpetaporNombre;
cabecera={'Fichero','Real'};
nhoja=NumFrames;
xlswrite(FicheroResumen,cabecera,nhoja);
xlswrite(FicheroResumen,NombreFichero',nhoja,'A2');
xlswrite(FicheroResumen,ClaseSonidoFichero',nhoja,'B2');
%}

%{
%-------------- Calcula la clase de sonido por ÁRBOL DE DECISIÓN
algoritmo=1;
cabecera={'Arbol Decision'};
display(cabecera);
clasificador=fitctree(matpar,ClasePatronNum);
[ClaseSonidoFichero,ProbClasFichero]=ClasificaCarpeta(algoritmo,clasificador);
EscribeColumnas(FicheroResumen,NumFrames,algoritmo,...
    cabecera,ClaseSonidoFichero,ProbClasFichero);
%}

%{
%-------------- Calcula la clase de sonido por kNN (k vecinos más próximos)
algoritmo=2;
cabecera={'kNN'};
display(cabecera);
clasificador=fitcknn(matpar,ClasePatronNum);
[ClaseSonidoFichero,ProbClasFichero]=ClasificaCarpeta(algoritmo,clasificador);
EscribeColumnas(FicheroResumen,NumFrames,algoritmo,...
    cabecera,ClaseSonidoFichero,ProbClasFichero);
%}

%{
%-------------- Calcula la clase de sonido por SVM (máquina de vectores soporte)
algoritmo=3;
cabecera={'SVM'};
display(cabecera);
clasificador=fitcMsvm(matpar,ClasePatronNum);
[ClaseSonidoFichero,ProbClasFichero]=ClasificaCarpeta(algoritmo,clasificador);
EscribeColumnas(FicheroResumen,NumFrames,algoritmo,...
    cabecera,ClaseSonidoFichero,ProbClasFichero);
%}

%{
%-------------- Calcula la clase de sonido por REGRESIÓN LOGÍSTICA
algoritmo=4;
cabecera={'Regresión logística'};
display(cabecera);
clasificador=mnrfit(matpar,ClasePatronNum);
[ClaseSonidoFichero,ProbClasFichero]=ClasificaCarpeta(algoritmo,clasificador);
EscribeColumnas(FicheroResumen,NumFrames,algoritmo,...
    cabecera,ClaseSonidoFichero,ProbClasFichero);
%}

%{
%-------------- Calcula la clase de sonido por redes neuronales
algoritmo=5;
cabecera={'Red neuronal'};
display(cabecera);
clasificador=feedforwardnet(10);
clasificador.trainParam.showWindow=0;
clasificador=train(clasificador,matpar',ClasePatronNum');
[ClaseSonidoFichero,ProbClasFichero]=ClasificaCarpeta(algoritmo,clasificador);
EscribeColumnas(FicheroResumen,NumFrames,algoritmo,...
    cabecera,ClaseSonidoFichero,ProbClasFichero);
%}

%{
%-------------- Calcula la clase de sonido por función discriminante
algoritmo=6;
cabecera={'Función discriminante'};
display(cabecera);
clasificador=fitcdiscr(matpar,ClasePatronNum);
[ClaseSonidoFichero,ProbClasFichero]=ClasificaCarpeta(algoritmo,clasificador);
EscribeColumnas(FicheroResumen,NumFrames,algoritmo,...
    cabecera,ClaseSonidoFichero,ProbClasFichero);
%}

%{
%-------------- Calcula la clase de sonido por redes bayesianas
algoritmo=7;
cabecera={'Red Bayesiana'};
display(cabecera);
clasificador=fitNaiveBayes(matpar,ClasePatronNum);
[ClaseSonidoFichero,ProbClasFichero]=ClasificaCarpeta(algoritmo,clasificador);
EscribeColumnas(FicheroResumen,NumFrames,algoritmo,...
    cabecera,ClaseSonidoFichero,ProbClasFichero);
%}


%{
%-------------- Calcula la clase de sonido por distancia mínima
algoritmo=8;
cabecera={'Distancia mínima'};
display(cabecera);
clasificador=fitDistMin(matpar,ClasePatronNum);
[ClaseSonidoFichero,ProbClasFichero]=ClasificaCarpeta(algoritmo,clasificador);
EscribeColumnas(FicheroResumen,NumFrames,algoritmo,...
    cabecera,ClaseSonidoFichero,ProbClasFichero);
%}

%{
%-------------- Calcula la clase de sonido por máxima verosimilitud
algoritmo=9;
cabecera={'Máxima Verosimilitud'};
display(cabecera);
clasificador=fitMaxVer(matpar,ClasePatronNum);
[ClaseSonidoFichero,ProbClasFichero]=ClasificaCarpeta(algoritmo,clasificador);
EscribeColumnas(FicheroResumen,NumFrames,algoritmo,...
    cabecera,ClaseSonidoFichero,ProbClasFichero);
%}


%
%-------------- Gráficos algoritmos
nhoja=NumFrames;
[~,~,hoja]=xlsread(FicheroResumen,nhoja);
tamhoja=size(hoja);
nfil=tamhoja(1);
ncol=tamhoja(2);
MatClas=cell2mat(hoja(2:nfil,2:ncol));
cabecera=hoja(1,:);
algoritmo=1;

GraficaPerfilExito(MatClas(:,1),MatClas(:,2*algoritmo))
%GraficaProbExito(MatClas(:,1),MatClas(:,2*algoritmo),MatClas(:,2*algoritmo+1))
%GraficaCompAlg(MatClas,cabecera);
%MatrizConfusion(hoja,algoritmo)
%EvaluacionClasificacion(hoja);
%GraficaROC(hoja);
%CohenKappaAnalysis(hoja);
%GraficaCompNumFrame(FicheroResumen);
%GraficaROCNumFrame(FicheroResumen);

kk=1;
%}


%
% Función que genera los datos de los patrones
%
function [X,Y,NomPar]=GeneraDatosPatrones
% X: matriz de predictores (observación,parámetro)
% Y: vector de clases (observación)
% NomPar: nombre de paráemtros (parámetro)

global PATRON;

[FicheroPatron,roiPatron,FicheroResumenPatrones]=DefinePatron;
nesp=PATRON.nesp;
X=[];
ClasePatronTxt=[];
for iesp=1:nesp
    FicheroMuestras=FicheroPatron{iesp};
    nmuest=length(FicheroMuestras);
    for imuest=1:nmuest
        NombreFichero=FicheroMuestras{imuest};
        [~,~,hoja]=xlsread(NombreFichero);
        tamhoja=size(hoja);
        nfil=tamhoja(1);
        ncol=tamhoja(2);
        matparx=cell2mat(hoja(2:nfil,2:ncol-1));
        ClasePatronTxtX=hoja(2:nfil,ncol);
        X=[X;matparx];
        ClasePatronTxt=[ClasePatronTxt;ClasePatronTxtX];
    end %de una muestra
end % de una especie
Y=ClaseTxt2Num(ClasePatronTxt);
npar=(ncol-2)/2;
X=X(:,1:npar);
NomPar=hoja(1,2:npar+1);


%
% Función que clasifica los sonidos de un fichero 
%   (a partir del fichero Excel)
%
function [ClaseSonido,ProbClase]=ClasificaFichero...
    (NombreFichero,algoritmo,clasificador)
% NombreFichero: nombre del fichero a procesar (sin extensión)
% algoritmo: Identificador del algoritmo de clasificación
% clasificador: Identificador del clasificador utilizado
% ClaseSonido: Resultado de la clasificación
% ProbClase: Probabilidad de la clasificación (factor de mérito)

global PATRON;
global ParUtil; %Parámetros utilizados
global NumPar NumFrames;

nesp=PATRON.nesp;
[~,~,hoja]=xlsread(NombreFichero);
tamhoja=size(hoja);
nfil=tamhoja(1);
ncol=tamhoja(2);
matparx=cell2mat(hoja(2:nfil,2:ncol-1));
nobservaciones=size(matparx,1);
matpar=zeros(nobservaciones-NumFrames+1,NumFrames*NumPar);
for k=NumFrames:nobservaciones
    matx=[];
    k0=k-NumFrames+1;
    for k1=k0:k
        x=matparx(k1,ParUtil);
        matx=[matx,x];
    end
    matpar(k0,:)=matx;
end

switch algoritmo
    case 1 % Árbol de decisión
        [PredSonido,ProbSonido]=predict(clasificador,matpar);
        
    case 2 % kNN
        [PredSonido,ProbSonido]=predict(clasificador,matpar);

    case 3 % SVM
        [PredSonido,ProbSonido]=svmPredict(clasificador,matpar);

    case 4 % Regresión logística
        ProbSonido= mnrval(clasificador,matpar);
        [~,PredSonido]=max(ProbSonido,[],2);

    case 5 % Red neuronal
        PredSonido=round(clasificador(matpar')');
        PredSonido=min(PredSonido,nesp+1);
        PredSonido=max(PredSonido,1);
        ProbSonido=zeros(nobservaciones,nesp+1);
        for k=1:nobservaciones-NumFrames+1
            ProbSonido(k,PredSonido(k))=1;
        end
        
    case 6 % Función discriminante
        [PredSonido,ProbSonido]=predict(clasificador,matpar);

    case 7 % Red bayesiana
        [ProbSonido,PredSonido]=posterior(clasificador,matpar);

    case 8 % Distancia mínima
        [PredSonido,ProbSonido]=DistMinPredict(clasificador,matpar);

    case 9 % Máxima versoimilitud
        [PredSonido,ProbSonido]=MaxVerPredict(clasificador,matpar);
end

ContEsp=zeros(1,nesp);
for iesp=1:nesp
    ResComp=eq(PredSonido,iesp*ones(length(PredSonido),1));
    ContEsp(iesp)=sum(ResComp);
end
[~,ClaseSonido]=max(ContEsp);
ProbClase=ContEsp(ClaseSonido)/sum(ContEsp);
kk=1;


%
% Función que clasifica los sonidos de una carpeta 
%   (a partir de los nombres de los ficheros Excel)
%
function [NombreFichero,ClaseSonidoFichero]=ClasificaCarpetaporNombre
% NombreFichero: Nombre de cada fichero excel (vector)
% ClaseSonidoFichero: Clase de sonido de cada fichero excel (vector)

global carpeta; %Carpeta de trabajo de los sonidos

NombreFichero=ObtieneFicherosExcelCarpeta(carpeta);
nfich=length(NombreFichero);
for ifich=1:nfich
    nombre=NombreFichero{ifich};

    for i=length(nombre):-1:1
        car=nombre(i);
        if strcmp(car,'\')
            idSon=nombre(i+1:i+2);
            break;
        end
    end
    
    switch idSon
        case 'A1'
            ClaseSonidoFichero(ifich)=1;
        case 'A2'
            ClaseSonidoFichero(ifich)=1;
        case 'A3'
            ClaseSonidoFichero(ifich)=2;
        case 'B1'
            ClaseSonidoFichero(ifich)=3;
        case 'B2'
            ClaseSonidoFichero(ifich)=3;
    end

end %del fichero


%
% Función que clasifica los sonidos de una carpeta 
%   (a partir de los ficheros Excel)
%
function [ClaseSonidoFichero,ProbClasFichero]=...
    ClasificaCarpeta(algoritmo,clasificador)
% algoritmo: Identificador del algoritmo de clasificación
% clasificador: Identificador del clasificador utilizado
% ClaseSonidoFichero: Clase de sonido predicho de cada fichero excel (vector)
% ProbClasFichero: Probabilidad de la clase de sonido de cada fichero excel (vector)

global carpeta; %Carpeta de trabajo de los sonidos

NombreFichero=ObtieneFicherosExcelCarpeta(carpeta);
nfich=length(NombreFichero);
for ifich=1:nfich
    nombre=NombreFichero{ifich};
    disp(nombre);
    [ClaseSonidoFichero(ifich),ProbClasFichero(ifich)]=...
        ClasificaFichero(nombre,algoritmo,clasificador);
    disp(ClaseSonidoFichero(ifich));
end %del fichero


%
% Función que escribe las columnas del fichero excel correspondientes a un
%   algoritmo de clasificación
%
function EscribeColumnas(FicheroResumen,hoja,algoritmo,...
    cabecera,ClaseSonido,ProbClas)
% FicheroResumen: Nombre del fichero excel en el que se escribe
% hoja: Número de hoja en la que escribir
% algoritmo: Identificador del algoritmo de clasificación
% cabecera: Texo con la cabecera de las columnas
% ClaseSonido: Resultado de la clasificación
% ProbClas: Probabilidad de la clasificación (factor de mérito)

columna=2*algoritmo+1;
LetCol1=char('A'+columna-1);
LetCol2=char('A'+columna);
xlswrite(FicheroResumen,cabecera,hoja,[LetCol1,'1']);
xlswrite(FicheroResumen,ClaseSonido',hoja,[LetCol1,'2']);
xlswrite(FicheroResumen,ProbClas',hoja,[LetCol2,'2']);

%
% Función que traza el gráfico de perfíl de éxito de una técnica de
% clasificación
%
function GraficaPerfilExito(SonidoReal,ClaseSonido)
% ClaseSonido: Clase de sonido predicho de cada fichero excel (vector)
% SonidoReal: Clase de sonido real de cada fichero excel (vector)

global PATRON;

nesp=PATRON.nesp;
nfichesp=zeros(nesp,1);
nfich=length(ClaseSonido);
for ifich=1:nfich
    nfichesp(SonidoReal(ifich))=nfichesp(SonidoReal(ifich))+1;
end

confusion=confusionmat(SonidoReal,ClaseSonido);
P=sum(confusion,2); %Positive
TP=diag(confusion); %True Positive
FN=P-TP; %False Negative

PP=sum(confusion,1)'; %Predicted Positive
FP=PP-TP; %False Positive
N=sum(P)-P; %Negative
TN=N-FP; %True Negative

% Accuracy
ACCx=(TP+TN)./(TP+TN+FP+FN); %Accuracy para cada clase
ACC = mean(ACCx); % Accuracy para el conjunto de las clases

% Gráfico del perfil de éxito
theta=linspace(0,2*pi,nesp+1);
tasaexito=zeros(1,nesp+1);
tasaexito(1:nesp)=TP./nfichesp;
tasaexito(nesp+1)=tasaexito(1);
rectangle('Position',[-1,-1,2,2],'Curvature',[1,1],'FaceColor','w');
axis off;
axis equal;
for tasa=[25,50,75]/100
    rectangle('Position',[-tasa,-tasa,2*tasa,2*tasa],'Curvature',[1,1],...
        'LineStyle',':');
end
for iesp=1:nesp
    th1=theta(iesp);
    xf=cos(th1);
    yf=sin(th1);
    color=PATRON.coloresp(iesp,:);
    line([0,xf],[0,yf],'Color',color);
    
    x1=tasaexito(iesp)*cos(th1);
    y1=tasaexito(iesp)*sin(th1);
    th2=theta(iesp+1);
    x2=tasaexito(iesp+1)*cos(th2);
    y2=tasaexito(iesp+1)*sin(th2);
    line([x1,x2],[y1,y2],'Color','k');
    
    delta=0.1;
    x1=x1-delta/2;
    y1=y1-delta/2;
    rectangle('Position',[x1,y1,delta,delta],'Curvature',[1,1],...
        'FaceColor',color);

    if xf>0
        ha='Left';
    elseif xf<0
        ha='Right';
    else
        ha='Center';
    end
    rt=1.1;
    xt=rt*cos(th1);
    yt=rt*sin(th1);
    label=[PATRON.NombreEspecie{iesp},' (',num2str(100*tasaexito(iesp),'%3.0f'),'%)'];
    text(xt,yt,label,'HorizontalAlignment',ha,'Color',color);
end
axis([-2,2,-1,1]);

exitoglobal=sum(TP)/nfich*100;
%title(['Exito global: ',num2str(exitoglobal,'%6.2f'),'%']);
%title(['Overall success rate: ',num2str(exitoglobal,'%6.2f'),'%']);
title(['Accuracy: ',num2str(100*ACC,'%6.2f'),'%']);




%{
%
% Función que traza el gráfico de perfíl de éxito de una técnica de
% clasificación
%
function GraficaPerfilExito(SonidoReal,ClaseSonido)
% ClaseSonido: Clase de sonido predicho de cada fichero excel (vector)
% SonidoReal: Clase de sonido real de cada fichero excel (vector)

global PATRON;

nesp=PATRON.nesp;
nexito=zeros(1,nesp);
nfichesp=zeros(1,nesp);
nfich=length(ClaseSonido);
for ifich=1:nfich
    nfichesp(SonidoReal(ifich))=nfichesp(SonidoReal(ifich))+1;
    if ClaseSonido(ifich)==SonidoReal(ifich)
        nexito(ClaseSonido(ifich))=nexito(ClaseSonido(ifich))+1;
    end
end

% Gráfico del perfil de éxito
theta=linspace(0,2*pi,nesp+1);
tasaexito=zeros(1,nesp+1);
tasaexito(1:nesp)=nexito./nfichesp;
tasaexito(nesp+1)=tasaexito(1);
rectangle('Position',[-1,-1,2,2],'Curvature',[1,1],'FaceColor','w');
axis off;
axis equal;
for tasa=[25,50,75]/100
    rectangle('Position',[-tasa,-tasa,2*tasa,2*tasa],'Curvature',[1,1],...
        'LineStyle',':');
end
for iesp=1:nesp
    th1=theta(iesp);
    xf=cos(th1);
    yf=sin(th1);
    color=PATRON.coloresp(iesp,:);
    line([0,xf],[0,yf],'Color',color);
    
    x1=tasaexito(iesp)*cos(th1);
    y1=tasaexito(iesp)*sin(th1);
    th2=theta(iesp+1);
    x2=tasaexito(iesp+1)*cos(th2);
    y2=tasaexito(iesp+1)*sin(th2);
    line([x1,x2],[y1,y2],'Color','k');
    
    delta=0.1;
    x1=x1-delta/2;
    y1=y1-delta/2;
    rectangle('Position',[x1,y1,delta,delta],'Curvature',[1,1],...
        'FaceColor',color);

    if xf>0
        ha='Left';
    elseif xf<0
        ha='Right';
    else
        ha='Center';
    end
    rt=1.1;
    xt=rt*cos(th1);
    yt=rt*sin(th1);
    label=[PATRON.NombreEspecie{iesp},' (',num2str(100*tasaexito(iesp),'%3.0f'),'%)'];
    text(xt,yt,label,'HorizontalAlignment',ha,'Color',color);
end
axis([-2,2,-1,1]);

exitoglobal=sum(nexito)/nfich*100;
%title(['Exito global: ',num2str(exitoglobal,'%6.2f'),'%']);
title(['Overall success rate: ',num2str(exitoglobal,'%6.2f'),'%']);
%}


%
% Función que traza el gráfico de la probabilidad de éxito 
%   de una técnica de clasificación
%
function GraficaProbExito(ClaseSonido,SonidoReal,ProbClase)
% ClaseSonido: Clase de sonido predicho de cada fichero excel (vector)
% SonidoReal: Clase de sonido real de cada fichero excel (vector)
% ProbClase: Probabilidad de la clase de sonido real de cada fichero excel (vector)

global PATRON;

nesp=PATRON.nesp;
nfichesp=zeros(1,nesp);
nfich=length(SonidoReal);
for ifich=1:nfich
    nfichesp(ClaseSonido(ifich))=nfichesp(ClaseSonido(ifich))+1;
end

% Dibujo de los áreas de cada especie real
xi=0;
for iesp=1:nesp
    if iesp~=nesp
        xf=sum(nfichesp(1:iesp))+0.5;
    else
        xf=sum(nfichesp(1:iesp))+1;
    end
    color=PATRON.coloresp(iesp,:);
    patch([xi,xf,xf,xi],[0,0,100,100],color,...
        'EdgeColor','none','FaceColor',color,'FaceAlpha',0.2);
    hold on;
    xi=xf;
end
axis([0,nfich+1,0,100]);

% Dibujo de la decisión de cada fichero
for ifich=1:nfich
    color=PATRON.coloresp(SonidoReal(ifich),:);
    plot([ifich,ifich],[0,100*ProbClase(ifich)],'Color',color);
    hold on;
    plot([ifich,ifich],100*[ProbClase(ifich),ProbClase(ifich)],'.','Color',color);
end %del fichero


%
% Función que traza el gráfico de comparación de número de Frames
%
function GraficaCompNumFrame(FicheroResumen)
% FicheroResumen: Nombre del fichero con el resumen de cálculos

global PATRON;
global ParUtil; %Parámetros utilizados
global NumPar NumFrames;

nesp=PATRON.nesp;
[~,~,hoja]=xlsread(FicheroResumen,1);
nfil=size(hoja,1);
ncol=size(hoja,2);
MatClas=cell2mat(hoja(2:nfil,2:ncol));
nfich=nfil-1;
ClaseReal=MatClas(:,1);
nomAlg={'ArDec','kNN','SVM','RegLog','Neur','Discr','Bayes','DisMin','MaxVer'};
nalg=length(nomAlg);
%nalg=(ncol-2)/2;
ordalg=[8,9,1:7];
Distancia=zeros(NumFrames,nalg);
Error=zeros(NumFrames,nalg);
for inumframes=1:NumFrames
    [~,~,hoja]=xlsread(FicheroResumen,inumframes);
    ncol=size(hoja,2);
    MatClas=cell2mat(hoja(2:nfil,2:ncol));
    nfichesp=zeros(nesp,1);
    nexito=zeros(nesp,nalg);
    %Cálculo del número de aciertos
    for ifich=1:nfich
        nfichesp(ClaseReal(ifich))=nfichesp(ClaseReal(ifich))+1;
        for ialg=1:nalg
            colalg=2*ialg;
            ClasePred=MatClas(:,colalg);
            if ClasePred(ifich)==ClaseReal(ifich)
                nexito(ClaseReal(ifich),ialg)=nexito(ClaseReal(ifich),ialg)+1;
            end
        end
    end
    %Cálculo de las tasas de éxito
    TasaExito=zeros(nesp,nalg);
    TasaExitoGlobal=zeros(1,nalg);
    RangoExitoGlobal=zeros(1,nalg);
    for ialg=1:nalg
        TasaExito(:,ialg)=100*nexito(:,ialg)./nfichesp;
        TasaExitoGlobal(ialg)=100*sum(nexito(:,ialg))./nfich;
        RangoExitoGlobal(ialg)=max(TasaExito(:,ialg))-min(TasaExito(:,ialg));
        x=100-TasaExitoGlobal(ialg);
        y=RangoExitoGlobal(ialg);
        Distancia(inumframes,ialg)=sqrt(x^2+y^2);
        Error(inumframes,ialg)=x;
    end
end %del número de parámetros
Merito=100-Distancia/sqrt(2);

%
% Gráfico de la tasa de éxito de cada algoritmo
MarkerVec=['+','o','*','s','d','p'];
xText=cell(nalg,1);
for ix=1:nalg
    ialg=ordalg(ix);
    xText(ix)=nomAlg(ialg);
    imarker=rem(ix-1,length(MarkerVec))+1;
    plot(1:NumFrames,Merito(:,ialg),...
        'Marker',MarkerVec(imarker));
    hold all;
end
%xlabel('Tamaño de la ventana');
%ylabel('Mérito (%)');
xlabel('Window size');
ylabel('Figure of merit (%)');
legend(xText,'Location','BestOutside');

figure;
maxMerito=max(Merito,[],2);
meanMerito=mean(Merito,2);
minMerito=min(Merito,[],2);
plot(1:NumFrames,meanMerito);
hold all;
plot(1:NumFrames,minMerito);
plot(1:NumFrames,maxMerito);
%xlabel('Tamaño de la ventana');
%ylabel('Mérito (%)');
xlabel('Window size');
ylabel('Figure of merit (%)');
ylim([0,120]);
%leyenda=[{'Máximo'},{'Media'},{'Mínimo'}];
leyenda=[{'Mean'},{'Min'},{'Max'}];
legend(leyenda,'Location','NorthEast');

figure;
dmaxMerito=max(Merito,[],2)-max(Merito(1,:));
plot(1:NumFrames,dmaxMerito);
%xlabel('Tamaño de la ventana');
%ylabel('\Delta Mérito (%)');
xlabel('Window size');
ylabel('\Delta figure of merit (%)');

figure;
dminError=min(Error,[],2)-min(Error(1,:));
plot(1:NumFrames,dminError);
%xlabel('Tamaño de la ventana');
%ylabel('\Delta Mérito (%)');
xlabel('Window size');
ylabel('\Delta error rate (%)');


%
% Función que traza el gráfico de comparación del ROC vs. número de Frames
%
function GraficaROCNumFrame(FicheroResumen)
% FicheroResumen: Nombre del fichero con el resumen de cálculos

global PATRON;
global ParUtil; %Parámetros utilizados
global NumPar NumFrames;

nomAlg={'DisMin','MaxVer','ArDec','kNN','SVM','RegLog','Neur','Discr','Bayes'};
nalg=length(nomAlg);
ROC=zeros(NumFrames,nalg);
F1=zeros(NumFrames,nalg);
%Cálculo de la evaluación ROC
for inumframes=1:NumFrames
    [~,~,hoja]=xlsread(FicheroResumen,inumframes);
    performance=EvaluacionClasificacion(hoja);
    %TPR=performance(:,4);
    %FPR=1-performance(:,5);
    %ROC(inumframes,:)=100*sqrt( FPR.^2 + (1-TPR).^2 ) / sqrt(2);
    PRC=performance(:,3);
    SNS=performance(:,4);
    SPC=1-performance(:,5);
    ROC(inumframes,:)=100*sqrt( SNS.^2 + SPC.^2 ) / sqrt(2);
    F1(inumframes,:)= 200* (PRC.*SNS) ./ (PRC+SNS); %F1 score
end %del número de parámetros

%
% Gráfico de la tasa de éxito de cada algoritmo
MarkerVec=['+','o','*','s','d','p'];
xText=cell(nalg,1);
%par = ROC;
par = F1;
for ialg=1:nalg
    xText(ialg)=nomAlg(ialg);
    imarker=rem(ialg-1,length(MarkerVec))+1;
    plot(1:NumFrames,par(:,ialg),...
        'Marker',MarkerVec(imarker));
    hold all;
end
%xlabel('Tamaño de la ventana');
xlabel('Window size');
%ylabel('ROC (%)');
ylabel('F_1 (%)');
legend(xText,'Location','BestOutside');

figure;
%par = ROC;
par = F1;
plot(1:NumFrames,max(par,[],2));
hold all;
plot(1:NumFrames,mean(par,2));
algopt=3; %Algoritmo óptimo 
%plot(1:NumFrames,par(:,algopt));
plot(1:NumFrames,min(par,[],2));
xlabel('Window size');
%ylabel('ROC (%)');
ylabel('F_1 (%)');
%leyenda=[{'Max'},{'Mean'},{'Best'},{'Min'}];
leyenda=[{'Max'},{'Mean'},{'Min'}];
legend(leyenda,'Location','SouthEast');
ylim([0,100]);

figure;
algopt=3;
%dminROC=min(ROC,[],2)-min(ROC(1,:));
dminROC=ROC(:,algopt)-ROC(1,algopt);
plot(1:NumFrames,dminROC);
%xlabel('Tamaño de la ventana');
%ylabel('\Delta ROC (%)');
xlabel('Window size');
ylabel('\Delta ROC (%)');


%
% Función que calcula la matriz de confusion de una técnica de clasificación
%
function MatrizConfusion(hoja,algoritmo)
% hoja: Hoja con datos de clase de sonido para cada fichero excel y cada algoritmo
% algoritmo: identificador del algoritmo del que se obtiene la matriz de confusión

global PATRON;
tamhoja=size(hoja);
nfil=tamhoja(1);
ncol=tamhoja(2);
MatClas=cell2mat(hoja(2:nfil,2:ncol));
SonidoReal=MatClas(:,1);
ClaseSonido=MatClas(:,2*algoritmo);

nesp=PATRON.nesp;

% Cálculo de la matriz de confusión (valores absolutos)
confusion=confusionmat(SonidoReal,ClaseSonido);
nfichesp=sum(confusion,2);
display(confusion);

% Cálculo de la matriz de confusión (porcentajes)
for iesp=1:nesp
    confusion(iesp,:)=confusion(iesp,:) / nfichesp(iesp)*100;
end
display(confusion);


%
% Función que calcula la evaluación de las técnicas de clasificación
%
function performance=EvaluacionClasificacion(hoja)
% hoja: Hoja con datos de clase de sonido para cada fichero excel y cada algoritmo
% performance: matriz de prestaciones (algoritmo,indicador)
%               1: Accuracy
%               2: Error Rate
%               3: Precision
%               4: Sensitivity
%               5: Specificity
%               6: ROC
%               7: F1 score

tamhoja=size(hoja);
nfil=tamhoja(1);
ncol=tamhoja(2);
MatClas=cell2mat(hoja(2:nfil,2:ncol));
SonidoReal=MatClas(:,1);

%nomAlg={'ArDec','kNN','SVM','RegLog','Neur','Discr','Bayes','DisMin','MaxVer'};
nomAlg={'DecTr','kNN','SVM','LogReg','Neur','Discr','Bayes','MinDis','MaxLik'};
ordalg=[8,9,1:7];
nalg=length(ordalg);
performance=zeros(nalg,5);
RowText = cell(1,nalg);
for ix=1:nalg
    ialg=ordalg(ix);
    ClaseSonido=MatClas(:,2*ialg);
    RowText{ix}=nomAlg{ialg};

    % Cálculo de la matriz de confusión (valores absolutos)
    confusion=confusionmat(SonidoReal,ClaseSonido);
    P=sum(confusion,2); %Positive
    TP=diag(confusion); %True Positive
    FN=P-TP; %False Negative

    PP=sum(confusion,1)'; %Predicted Positive
    FP=PP-TP; %False Positive
    N=sum(P)-P; %Negative
    TN=N-FP; %True Negative

    % Prestaciones para cada clase
    ACCx=(TP+TN)./(TP+TN+FP+FN); %Accuracy
    ERRx=(FP+FN)./(TP+TN+FP+FN); %Error Rate
    PRCx=TP./(TP+FP); %Precision
    SNSx=TP./(TP+FN); %Sensitivity
    SPCx=TN./(FP+TN); %Specificity
    PRCx(isnan(PRCx)) = 0; %Elimina NaN

    % Prestaciones para el conjunto de las clases
    ACCm = mean(ACCx);
    ERRm = mean(ERRx);
    PRCm = mean(PRCx);
    SNSm = mean(SNSx);
    SPCm = mean(SPCx);
    ROCm = sqrt( SNSm.^2 + SPCm.^2 ) / sqrt(2); %ROC
    F1m = 2* (PRCm.*SNSm) ./ (PRCm+SNSm); %F1 score

    performance(ix,1)=ACCm; %Accuracy
    performance(ix,2)=ERRm; %Error Rate
    performance(ix,3)=PRCm; %Precision
    performance(ix,4)=SNSm; %Sensitivity
    performance(ix,5)=SPCm; %Specificity
    performance(ix,6)=ROCm; %ROC
    performance(ix,7)=F1m; %F1 score
end
ACC = performance(:,1);
ERR = performance(:,2);
PRC = performance(:,3);
SNS = performance(:,4);
SPC = performance(:,5);
ROC = performance(:,6);
F1 = performance(:,7);
T = table(ACC,ERR,PRC,SNS,SPC,ROC,F1,'RowNames',RowText);
disp(T);
kk=1;


%
% Función que traza el gráfico de comparación de algoritmos
%
function GraficaCompAlg(hoja,cabecera)
% hoja: Hoja con datos de clase de sonido para cada fichero excel y cada algoritmo
% cabecera: Cabecera de las columnas de las hojas (vector)

global PATRON;

nesp=PATRON.nesp;
tamhoja=size(hoja);
nfil=tamhoja(1);
ncol=tamhoja(2);
nfich=nfil;
nomAlg={'ArDec','kNN','SVM','RegLog','Neur','Discr','Bayes','DisMin','MaxVer'};
ordalg=[8,9,1:7];
nalg=length(ordalg);
ClaseReal=hoja(:,1);
nfichesp=zeros(nesp,1);
nexito=zeros(nesp,nalg);
for ifich=1:nfich
    nfichesp(ClaseReal(ifich))=nfichesp(ClaseReal(ifich))+1;
    for ix=1:nalg
        ialg=ordalg(ix);
        colalg=2*ialg;
        ClasePred=hoja(:,colalg);
        if ClasePred(ifich)==ClaseReal(ifich)
            nexito(ClaseReal(ifich),ix)=nexito(ClaseReal(ifich),ix)+1;
        end
    end
end

TasaExito=zeros(nesp,nalg);
TasaExitoGlobal=zeros(1,nalg);
RangoExitoGlobal=zeros(1,nalg);
xText=cell(nalg,1);
for ix=1:nalg
    TasaExito(:,ix)=100*nexito(:,ix)./nfichesp;
    TasaExitoGlobal(ix)=100*sum(nexito(:,ix))./nfich;
    RangoExitoGlobal(ix)=max(TasaExito(:,ix))-min(TasaExito(:,ix));
    xText{ix}=nomAlg{ordalg(ix)};
end

%
% Gráfico de la tasa de éxito de cada algoritmo
bar(1:nalg,TasaExitoGlobal,'c');
hold all;
for iesp=1:nesp
    plot(1:nalg,TasaExito(iesp,:),'*');
end
axis([0,nalg+1,0,110]);
set(gca,'XTickLabel',xText);
ylabel('Tasa de éxito');
%}

% Gráfico de la tasa de éxito vs. Rango de éxito
figure;
MarkerVec=['+','o','*','s','d','p'];
for ialg=1:nalg
    imarker=rem(ialg-1,length(MarkerVec))+1;
    plot(100-TasaExitoGlobal(ialg),RangoExitoGlobal(ialg),...
        'Marker',MarkerVec(imarker),'LineStyle','none');
    hold all;
end
legend(xText,'Location','NorthEast');
xlabel('Tasa de error');
ylabel('Rango de tasa de error');
axis([0,100,0,100]);


%
% Función que traza el gráfico ROC de comparación de algoritmos
%   ROC: Receiver Operating Characteristic
%
function GraficaROC(hoja)
% hoja: Hoja con datos de clase de sonido para cada fichero excel y cada algoritmo

ordalg=[8,9,1:7];
nomAlg={'DecTr','kNN','SVM','LogReg','Neur','Discr','Bayes','MinDis','MaxLik'};
MarkerVec=['+','o','*','s','d','p'];
performance=EvaluacionClasificacion(hoja);
nalg=length(nomAlg);
xText=cell(nalg,1);
for ialg=1:nalg
    xText{ialg}=nomAlg{ordalg(ialg)};
    %xText{ialg}=nomAlg{ialg};
    imarker=rem(ialg-1,length(MarkerVec))+1;
    x=100*(1-performance(ialg,5)); %Fall-out (1-Specificity)=False Positive Rate
    y=100*performance(ialg,4); %Sensitivity=True Positive Rate
    plot(x,y,...
        'Marker',MarkerVec(imarker),'LineStyle','none');
    hold all;
end
plot([0,100],[0,100],':k');
legend(xText,'Location','NorthEast');
xlabel('False Positive Rate (%)');
ylabel('True Positive Rate (%)');
axis([0,100,0,100]);


%
% Función que realiza la comparación de algoritmos mediante el coeficiente kappa de Cohen
%
function CohenKappaAnalysis(hoja)
% hoja: Hoja con datos de clase de sonido para cada fichero excel y cada algoritmo

global PATRON;

nclases=PATRON.nesp;
tamhoja=size(hoja);
nfil=tamhoja(1);
ncol=tamhoja(2);
MatClas=cell2mat(hoja(2:nfil,2:ncol));
%nalg=(ncol-2)/2;

nomAlg={'DecTr','kNN','SVM','LogReg','Neur','Discr','Bayes','MinDis','MaxLik'};
ordalg=[8,9,1:7];
nalg=length(nomAlg);
NomAlgOrd=cell(1,nalg);
for i=1:nalg
    NomAlgOrd{i}=nomAlg{ordalg(i)};
end

% Calcula la matriz kappa de Cohen
KappaMat=zeros(nalg,nalg);
for ix=1:nalg
    ialg=ordalg(ix);
    ClaseAlgI=MatClas(:,2*ialg);
    for jx=ix:nalg
        jalg=ordalg(jx);
        ClaseAlgJ=MatClas(:,2*jalg);
        nfich=length(ClaseAlgI);
        X=zeros(nclases,nclases);
        for ifich=1:nfich
            i=ClaseAlgI(ifich);
            j=ClaseAlgJ(ifich);
            X(i,j)=X(i,j)+1;
        end %de ifich
        [po,pe,k,sek,ci,km,ratio,var,z,p]=kappa(X);
        KappaMat(ix,jx)=min(k,1);
        KappaMat(jx,ix)=min(k,1);
    end %de jalg
end %de ialg

% Traza el gráfico de la matriz kappa
colormap summer;
mapacolor=colormap;
ncolor=size(mapacolor,1);
for x=1:nalg
    for y=1:nalg
        z=max(0,KappaMat(x,y));
        color=mapacolor(round(z*(ncolor-1)+1),:);
        rectangle('Position',[x-0.5,y-0.5,1,1],'FaceColor',color);
        ztexto=[num2str(100*z,'%2.0f'),'%'];
        text(x,y,ztexto,'HorizontalAlignment','center');
    end
end
axis([0.5,nalg+0.5,0.5,nalg+0.5]);
%----- Etiquetas del eje horizontal girados 30º
set(gca,'XTick',1:nalg,'XTickLabel','');
hx = get(gca,'XLabel');
set(hx,'Units','data'); 
pos = get(hx,'Position'); 
t=zeros(1,nalg);
for i = 1:nalg
    t(i)=text(i,pos(2),NomAlgOrd{i}); 
end 
set(t,'Rotation',30,'HorizontalAlignment','right')
%----- Etiquetas del eje vertical
set(gca,'YTick',1:nalg+1,'YTickLabel',NomAlgOrd);
colorbar('YTick',1:ncolor/4:ncolor+1,'YTickLabel',...
    {'0%','25%','50%','75%','100%'});

% Traza el gráfico de cajas de la kappa media de cada algoritmo
figure;
Kappax=KappaMat;
A=diag(ones(1,nalg));
Kappax(A==1)=nan;
boxplot(Kappax,NomAlgOrd,'labelorientation','inline');
ylabel('Cohen''s \kappa rate');


%
% Obtiene los nombres de los ficheros Excel de una carpeta
%
function NombreFichero=ObtieneFicherosExcelCarpeta(carpeta)
% carpeta: nombre de la carpeta a analizar
% NombreFichero: nombre de los ficheros (vector de celdas)

dirant=cd;
cd(carpeta);
directorio=dir;
cd(dirant);

nficheros=length(directorio);
for i=1:nficheros
    nombres{i}=directorio(i).name;
end

[NombresOrd,iNombresOrd]=sort(nombres);
nfich=0;
for i=3:nficheros % Evita el '.' y el '..'
    j=iNombresOrd(i);
    NombreCompleto=directorio(j).name;
    longnom=length(NombreCompleto);
    idSon=NombreCompleto(1:3);
    switch idSon
        case 'A1-'
            idSonCorrecto=1;
        case 'A2-'
            idSonCorrecto=1;
        case 'A3-'
            idSonCorrecto=1;
        case 'B1-'
            idSonCorrecto=1;
        case 'B2-'
            idSonCorrecto=1;
        otherwise
            idSonCorrecto=0;
    end
    for k=longnom:-1:1
        if NombreCompleto(k)=='.'
            extension=NombreCompleto(k+1:longnom);
            break;
        end
    end
    if idSonCorrecto && strcmp(extension,'xls')
        nfich=nfich+1;
        NombreFichero{nfich}=[carpeta,NombreCompleto(1:k-1)];
    end
end


%
% Función que define las variables generales del programa
%
function DefineVariablesGenerales

global GENERAL; % Estructura de datos generales
    %GENERAL.hopSize: Duración de un "salto" de tiempo en el espectrograma (MP7)
    %GENERAL.msegSize: Duración de un microsegmento
    %GENERAL.frPrelev: Frecuencias de la banda relevante de potencias [baja,alta]
    %GENERAL.nForm: Número máximo de formantes considerados
    %GENERAL.NomParFrame: Nombre de los parámetros de un frame
    %GENERAL.NomCorParFrame: Nombre corte de los parámetros de un frame
    %GENERAL.ParSigEsp: Nombre de los parámetros significativos para determinar una especie en un frame
    %GENERAL.NumParPrimarios: Número de parámetros de frame primarios
    %GENERAL.NumParSecundarios: Número de parámetros de frame secundarios
    %GENERAL.NumParTerciarios: Número de parámetros de frame terciarios
    %GENERAL.gradoLPC: grado del polinomio usado en el algoritmo LPC
    %GENERAL.frFilt: Frecuencias de corte del filtr de preprocesamiento[baja,alta]
global PATRON;
    

GENERAL.hopSize=0.01; %Duración de un "salto" de tiempo en el espectrograma (MP7)
GENERAL.msegSize=0.1; %Duración de un microsegmento
GENERAL.frPrelev=[500,5000]; %Frecuencias de la banda relevante de potencias
GENERAL.nForm=3; %Número de formantes
GENERAL.gradoLPC=30; %Grado del polinomio usado en el algoritmo LPC
GENERAL.frFilt=[300,10000]; %Frecuencias de corte del filtro de preprocesamiento[baja,alta]

nesp=PATRON.nesp;
nform=GENERAL.nForm;
i=0;
i=i+1; GENERAL.NomParFrame{i}='Pt';
i=i+1; GENERAL.NomParFrame{i}='PRelev';
i=i+1; GENERAL.NomParFrame{i}='CP';
i=i+1; GENERAL.NomParFrame{i}='DE';
i=i+1; GENERAL.NomParFrame{i}='Pl';
i=i+1; GENERAL.NomParFrame{i}='Pitch';
i=i+1; GENERAL.NomParFrame{i}='Ra';
i=i+1; GENERAL.NomParFrame{i}='Fla';
for k=1:nform
    i=i+1; GENERAL.NomParFrame{i}=['FreqForm',num2str(k)];
    i=i+1; GENERAL.NomParFrame{i}=['AbForm',num2str(k)];
end
i=i+1; GENERAL.NomParFrame{i}='CA';
i=i+1; GENERAL.NomParFrame{i}='DeA';
i=i+1; GENERAL.NomParFrame{i}='DiA';
i=i+1; GENERAL.NomParFrame{i}='VA';
GENERAL.NumParPrimarios=i;

i=i+1; GENERAL.NomParFrame{i}='PtSigma';
i=i+1; GENERAL.NomParFrame{i}='PRelevSigma';
i=i+1; GENERAL.NomParFrame{i}='CPSigma';
i=i+1; GENERAL.NomParFrame{i}='DESigma';
i=i+1; GENERAL.NomParFrame{i}='PlSigma';
i=i+1; GENERAL.NomParFrame{i}='PitchSigma';
i=i+1; GENERAL.NomParFrame{i}='RaSigma';
i=i+1; GENERAL.NomParFrame{i}='FlaSigma';
for k=1:nform
    i=i+1; GENERAL.NomParFrame{i}=['FreqFormSigma',num2str(k)];
    i=i+1; GENERAL.NomParFrame{i}=['AbFormSigma',num2str(k)];
end
i=i+1; GENERAL.NomParFrame{i}='CASigma';
i=i+1; GENERAL.NomParFrame{i}='DeASigma';
i=i+1; GENERAL.NomParFrame{i}='DiASigma';
i=i+1; GENERAL.NomParFrame{i}='VASigma';
GENERAL.NumParSecundarios=i-GENERAL.NumParPrimarios;

for k=1:nesp
    i=i+1; GENERAL.NomParFrame{i}=['VrEsp',num2str(k)];
    i=i+1; GENERAL.NomParFrame{i}=['DisEsp',num2str(k)];
end
i=i+1; GENERAL.NomParFrame{i}='Relev';
GENERAL.NumParTerciarios=i-GENERAL.NumParPrimarios-GENERAL.NumParSecundarios;

i=0;
i=i+1; GENERAL.NomCorParFrame{i}='Pt';
i=i+1; GENERAL.NomCorParFrame{i}='PR';
i=i+1; GENERAL.NomCorParFrame{i}='CP';
i=i+1; GENERAL.NomCorParFrame{i}='DE';
i=i+1; GENERAL.NomCorParFrame{i}='Pl';
i=i+1; GENERAL.NomCorParFrame{i}='Pitch';
i=i+1; GENERAL.NomCorParFrame{i}='Ra';
i=i+1; GENERAL.NomCorParFrame{i}='Fla';
for k=1:nform
    i=i+1; GENERAL.NomCorParFrame{i}=['FrF',num2str(k)];
    i=i+1; GENERAL.NomCorParFrame{i}=['AbF',num2str(k)];
end
i=i+1; GENERAL.NomCorParFrame{i}='CA';
i=i+1; GENERAL.NomCorParFrame{i}='DeA';
i=i+1; GENERAL.NomCorParFrame{i}='DiA';
i=i+1; GENERAL.NomCorParFrame{i}='VA';

i=i+1; GENERAL.NomCorParFrame{i}='Pt~';
i=i+1; GENERAL.NomCorParFrame{i}='PR~';
i=i+1; GENERAL.NomCorParFrame{i}='CP~';
i=i+1; GENERAL.NomCorParFrame{i}='DE~';
i=i+1; GENERAL.NomCorParFrame{i}='Pl~';
i=i+1; GENERAL.NomCorParFrame{i}='Pitch~';
i=i+1; GENERAL.NomCorParFrame{i}='Ra~';
i=i+1; GENERAL.NomCorParFrame{i}='Fla~';
for k=1:nform
    i=i+1; GENERAL.NomCorParFrame{i}=['FF~',num2str(k)];
    i=i+1; GENERAL.NomCorParFrame{i}=['AB~',num2str(k)];
end
i=i+1; GENERAL.NomCorParFrame{i}='CA~';
i=i+1; GENERAL.NomCorParFrame{i}='DeA~';
i=i+1; GENERAL.NomCorParFrame{i}='DiA~';
i=i+1; GENERAL.NomCorParFrame{i}='VA~';

for k=1:nesp
    i=i+1; GENERAL.NomCorParFrame{i}=['VrE',num2str(k)];
    i=i+1; GENERAL.NomCorParFrame{i}=['DiE',num2str(k)];
end
i=i+1; GENERAL.NomCorParFrame{i}='Relev';


%
% Función que define los patrones
%
function [FicheroPatron,roiPatron,FicheroResumenPatrones]=DefinePatron
% FicheroPatron: Nombre de los ficheros Patrón (celda de mxn)
% roiPatron: Regiones de Interés (ROI) de los ficheros Patrón (celda de mxn)
% FicheroResumenPatrones: Nombre del fichero resumen de patrones (global PATRON)

global PATRON; % Estructura de datos de patrones
% Datos de cada patron de especie (se aproxima por una mezcla de N gaussianas)
    %PATRON.nesp: Número de especies en los patrones
    %PATRON.NombreEspecie: Nombre de las especies en los patrones
    %PATRON.coloresp: Color usado para las especies (matriz [tipo región,especie])
    %PATRON.ng: Número de gaussianas en la mezcla de funciones de distribución
    %PATRON.nt: Número de frames en los patrones (matriz [especie,muestra])
   

PATRON.nesp=3;
%PATRON.NombreEspecie{1}='Sapo corredor';
%PATRON.NombreEspecie{2}='Sapo corredor: suelta';
%PATRON.NombreEspecie{3}='Sapo partero común';
PATRON.NombreEspecie{1}='Calamita: mating';
PATRON.NombreEspecie{2}='Calamita: release';
PATRON.NombreEspecie{3}='Alytes obstetricans';
PATRON.NombreEspecieLatin{1}='Epidalea Calamita';
PATRON.NombreEspecieLatin{2}='Epidalea Calamita: suelta';
PATRON.NombreEspecieLatin{3}='Alytes obstetricans';
PATRON.ng=2;
tablaColores=[0.00,0.00,1.00;...
              0.00,0.50,0.00;...
              1.00,0.00,0.00;...
              0.75,0.00,0.75];
PATRON.coloresp=tablaColores(1:PATRON.nesp+1,:);
FicheroResumenPatrones='sonidos\Procesados\Patrones\IANTA-PATRON';
FicheroPatron=cell(1,PATRON.nesp);
FicheroPatron{1}={'sonidos\Procesados\Patrones\0501_0-10',...
                  'sonidos\Procesados\Patrones\0707_0-10'};
FicheroPatron{2}={'sonidos\Procesados\Patrones\0503'};
FicheroPatron{3}={'sonidos\Procesados\Patrones\1325',...
                  'sonidos\Procesados\Patrones\1330'};

roiPatron=cell(1,PATRON.nesp);
roiPatron{1}={[  21,  78;
                127, 231;
                264, 366;
                395, 489;
                523, 628;
                664, 759;
                798, 900;
                933,1035],...
              [  20,  44;
                 73, 127;
                150, 194;
                221, 270;
                294, 339;
                362, 411;
                435, 484;
                507, 545;
                569, 611;
                635, 679;
                704, 747;
                771, 821;
                845, 887;
                910, 950;
                976,1018]};

roiPatron{2}={[  42,  58;
                128, 136;
                377, 388;
                745, 753;
                795, 812;
               1019,1034;
               1345,1356;
               1403,1418;
               1463,1479;
               1570,1577;
               1632,1642;
               1858,1865;
               1915,1925;
               2241,2251;
               2303,2313;
               2366,2380;
               2518,2531;
               2574,2582;
               2810,2817;
               2858,2873]};
              
roiPatron{3}={[  73,  89;
                150, 164;
                497, 513;
                876, 892;
               1433,1450;
               1808,1825;
               2206,2222;
               2556,2572;
               2933,2948;
               3406,3422;
               3677,3699;
               3721,3737],...
              [ 182, 194;
                263, 277;
                696, 711;
               1226,1239;
               1737,1751;
               2374,2388;
               2919,2934;
               3436,3451;
               4001,4015;
               4487,4501;
               4990,5005]};


%
% Función que convierte de la clase de un sonido de texto a número
%
function ClaseNum=ClaseTxt2Num(ClaseTxt)
% ClaseTxt: Clase de sonido en formato texto (vector)
% ClaseNum: Clase de sonido en formato número (vector)

global PATRON;
nesp=PATRON.nesp;
nelem=length(ClaseTxt);
ClaseNum=(nesp+1)*ones(nelem,1);
for k=1:nelem
    for iesp=1:nesp
        if strcmp(ClaseTxt(k),PATRON.NombreEspecie{iesp});
            ClaseNum(k)=iesp;
        end
    end
end


%
% Función que obtinen modelo SVM con múltiples clases
%
function clasificador=fitcMsvm(X,Y)
% X: matriz de predictores (observación,parámetro)
% Y: vector de clases (observación)
% clasificador: clasificador; vector (especies)

global PATRON;
nesp=PATRON.nesp;
nobservaciones=length(Y);
Y2=zeros(nobservaciones,1);
clasificador=cell(nesp+1,1);
for iesp=1:nesp+1
    for k=1:nobservaciones
        if Y(k)==iesp
            Y2(k)='A';
        else
            Y2(k)='B';
        end
    end %de las observaciones
    SVMModel=fitcsvm(X,Y2,'Standardize',true);
    CompactSVMModel=compact(SVMModel);
    clasificador{iesp}=fitPosterior(CompactSVMModel,X,Y2);
end %de las especies
save('SVM','clasificador');
kk=1;


%
% Función que clasifica mediante SVM con múltiples clases
%
function [Z,ProbZ]=svmPredict(clasificador,X)
% X: matriz de predictores (observación,parámetro)
% clasificador: clasificador; vector (clases)
% Z: vector de predicciones (observación)
% ProbZ: Probabilidad de las predicciones: matriz (observación,especie)

global PATRON;
nesp=PATRON.nesp;
nobservaciones=length(X);
ProbZ2=zeros(nobservaciones,2,nesp+1);
for iesp=1:nesp+1
    modeloSVM=clasificador{iesp};
    [~,ProbZ2(:,:,iesp)]=predict(modeloSVM,X);
end

ProbZ(:,:)=ProbZ2(:,1,:);
Z=zeros(nobservaciones,1);
for k=1:nobservaciones
    ProbZ(k,:)=ProbZ(k,:)/sum(ProbZ(k,:));
    [~,Z(k)]=max(ProbZ(k,:));
end
kk=1;


%
% Función que obtiene el clasificador por máxima verosimilitud
%
function clasificador=fitMaxVer(X,Y)
% X: matriz de predictores (observación,parámetro)
% Y: vector de clases (observación)
% clasificador: clasificador; vector (clases)

global PATRON;
nesp=PATRON.nesp;
ng=PATRON.ng;
npar=size(X,2);
gmd=cell(nesp+1);
for iesp=1:nesp+1
    Xesp=X(Y==iesp,:);
    GMDini.PComponents=0.5*ones(1,ng);
    medini=mean(Xesp);
    varini=var(Xesp);
    for i=1:npar
        GMDini.mu(:,i)=linspace(medini(i)-sqrt(varini(i)),medini(i)+sqrt(varini(i)),ng)';
    end
    for i=1:ng
        GMDini.Sigma(:,:,i)=cov(Xesp);
    end
    warning off all;
    %gmd{iesp}=fitgmdist(param,ng,'Start',GMDini);
    try
        gmd{iesp}=fitgmdist(Xesp,ng,'Start',GMDini);
    catch exception
        gmd{iesp}=gmdistribution(GMDini.mu,GMDini.Sigma,GMDini.PComponents);
    end
    warning on all;

end
clasificador=gmd;


%
% Función que clasifica mediante Máxima Verosimilitud
%
function [Z,ProbZ]=MaxVerPredict(clasificador,X)
% X: matriz de predictores (observación,parámetro)
% clasificador: clasificador
% Z: vector de predicciones (observación)
% ProbZ: Probabilidad de las predicciones: matriz (observación,especie)

global PATRON;
nesp=PATRON.nesp;
ng=PATRON.ng;
nobservaciones=size(X,1);
gmdEsp=clasificador;
Z=zeros(nobservaciones,1);
ProbZ=zeros(nobservaciones,nesp+1);
verosimilitud=zeros(nesp+1,1);

for k=1:nobservaciones
    fdp=zeros(1,nesp+1);
    x2=X(k,:);
    for iesp=1:nesp+1
        for in=1:ng
            gmd=gmdEsp{iesp};
            fdp(iesp)=fdp(iesp)+gmd.PComponents(in)*...
                mvnpdf(x2,gmd.mu(in,:),gmd.Sigma(:,:,in));
        end %de gaussiana
    verosimilitud(iesp)=log10(max(fdp(iesp),realmin));
    end %de especie
    [~,espopt]=max(verosimilitud);
    Z(k)=espopt;
    ProbZ(k,:)=verosimilitud;
end %de frame
kk=1;


%
% Función que obtiene el clasificador por distancia mínima
%
function clasificador=fitDistMin(X,Y)
% X: matriz de predictores (observación,parámetro)
% Y: vector de clases (observación)
% clasificador: clasificador; vector (clases)

global PATRON;
nesp=PATRON.nesp;
npar=size(X,2);
mu=zeros(nesp+1,npar);
sigma=zeros(nesp+1,npar);
for iesp=1:nesp+1
    Xesp=X(Y==iesp,:);
    mu(iesp,:)=median(Xesp,1);
    sigma(iesp,:)=iqr(Xesp,1);
end
clasificador={mu,sigma};


%
% Función que clasifica mediante Distancia Mínima
%
function [Z,ProbZ]=DistMinPredict(clasificador,X)
% X: matriz de predictores (observación,parámetro)
% clasificador: clasificador
% Z: vector de predicciones (observación)
% ProbZ: Probabilidad de las predicciones: matriz (observación,especie)

global PATRON;
nesp=PATRON.nesp;
nobservaciones=size(X,1);
mu=clasificador{1};
Z=zeros(nobservaciones,1);
ProbZ=zeros(nobservaciones,nesp+1);
distancia=zeros(nesp+1,1);
rangopar=max(X,[],1)-min(X,[],1);

for k=1:nobservaciones
    x2=X(k,:);
    for iesp=1:nesp+1
        mux=mu(iesp,:);
        distpar=abs(x2-mux);
        distnorm=distpar./rangopar;
        distancia(iesp)=sqrt(sum(distnorm.^2));
    end %de especie
    
    score=sum(distancia)./distancia;
    ProbZ(k,:)=score./sum(score);
    [~,espmin]=min(distancia);
    Z(k)=espmin;
end %de frame
kk=1;


function [po,pe,k,sek,ci,km,ratio,var,z,p]=kappa(varargin)
% KAPPA: This function computes the Cohen's kappa coefficient.
% Cohen's kappa coefficient is a statistical measure of inter-rater
% reliability. It is generally thought to be a more robust measure than
% simple percent agreement calculation since k takes into account the
% agreement occurring by chance.
% Kappa provides a measure of the degree to which two judges, A and B,
% concur in their respective sortings of N items into k mutually exclusive
% categories. A 'judge' in this context can be an individual human being, a
% set of individuals who sort the N items collectively, or some non-human
% agency, such as a computer program or diagnostic test, that performs a
% sorting on the basis of specified criteria.
% The original and simplest version of kappa is the unweighted kappa
% coefficient introduced by J. Cohen in 1960. When the categories are
% merely nominal, Cohen's simple unweighted coefficient is the only form of
% kappa that can meaningfully be used. If the categories are ordinal and if
% it is the case that category 2 represents more of something than category
% 1, that category 3 represents more of that same something than category
% 2, and so on, then it is potentially meaningful to take this into
% account, weighting each cell of the matrix in accordance with how near it
% is to the cell in that row that includes the absolutely concordant items.
% This function can compute a linear weights or a quadratic weights.
%
% Syntax: [po,pe,k,sek,ci,km,ratio,var,z,p]=kappa(X,W,ALPHA)
%      
%     Inputs:
%           X - square data matrix
%           W - Weight (0 = unweighted; 1 = linear weighted; 2 = quadratic
%           weighted; -1 = display all. Default=0)
%           ALPHA - default=0.05.
%
%     Outputs:
%           po: Observed agreement
%           pe: Random agreement
%           k: Cohen's kappa
%               k<=0: Poor agreement
%               0.0<k<=0.2: Slight agreement
%               0.2<k<=0.4: Fair agreement
%               0.4<k<=0.6: Moderate agreement
%               0.6<k<=0.8 Substantial agreement
%               0.8<k<=1.0: Perfect agreement
%           sekf: kappa error
%           ci: kappa confidence interval
%           km: Maximum possible kappa, given the observed marginal frequencies
%           ratio: kappa observed as proportion of maximum possible
%           var: variance
%           z: normalized kappa
%           p: observed agreement hypotesis test
%               p<0.05: Reject null hypotesis: observed agreement is not accidental
%               p>=0.05: Accept null hypotesis: observed agreement is accidental
%
%      Example: 
%
%           x=[88 14 18; 10 40 10; 2 6 12];
%
%           Calling on Matlab the function: kappa(x)
%
%           Answer is:
%
% UNWEIGHTED COHEN'S KAPPA
% --------------------------------------------------------------------------------
% Observed agreement (po) = 0.7000
% Random agreement (pe) = 0.4100
% Agreement due to true concordance (po-pe) = 0.2900
% Residual not random agreement (1-pe) = 0.5900
% Cohen's kappa = 0.4915
% kappa error = 0.0549
% kappa C.I. (alpha = 0.0500) = 0.3839     0.5992
% Maximum possible kappa, given the observed marginal frequencies = 0.8305
% k observed as proportion of maximum possible = 0.5918
% Moderate agreement
% Variance = 0.0031     z (k/sqrt(var)) = 8.8347    p = 0.0000
% Reject null hypotesis: observed agreement is not accidental
%
%           Created by Giuseppe Cardillo
%           giuseppe.cardillo-edta@poste.it
%
% To cite this file, this would be an appropriate format:
% Cardillo G. (2007) Cohens kappa: compute the Cohen's kappa ratio on a 2x2 matrix.   
% http://www.mathworks.com/matlabcentral/fileexchange/15365



%Input Error handling

args=cell(varargin);
nu=numel(args);
if isempty(nu)
    error('Warning: Matrix of data is missed...')
elseif nu>3
    error('Warning: Max three input data are required')
end
default.values = {[],0,0.05};
default.values(1:nu) = args;
[x w alpha] = deal(default.values{:});

if isempty(x)
    error('Warning: X matrix is empty...')
end
if isvector(x)
    error('Warning: X must be a matrix not a vector')
end
if ~all(isfinite(x(:))) || ~all(isnumeric(x(:)))
    error('Warning: all X values must be numeric and finite')
end   
if ~isequal(x(:),round(x(:)))
    error('Warning: X data matrix values must be whole numbers')
end

m=size(x);
if ~isequal(m(1),m(2))
    error('Input matrix must be a square matrix')
end
if nu>1 %eventually check weight
    if ~isscalar(w) || ~isfinite(w) || ~isnumeric(w) || isempty(w)
        error('Warning: it is required a scalar, numeric and finite Weight value.')
    end
    a=-1:1:2;
    if isempty(a(a==w))%check if w is -1 0 1 or 2
        error('Warning: Weight must be -1 0 1 or 2.')
    end
end
if nu>2 %eventually check alpha
    if ~isscalar(alpha) || ~isnumeric(alpha) || ~isfinite(alpha) || isempty(alpha)
        error('Warning: it is required a numeric, finite and scalar ALPHA value.');
    end
    if alpha <= 0 || alpha >= 1 %check if alpha is between 0 and 1
        error('Warning: ALPHA must be comprised between 0 and 1.')
    end
end
clear args default nu

m(2)=[]; 
if w==0 || w==-1
    f=diag(ones(1,m)); %unweighted
    [po,pe,k,sek,ci,km,ratio,var,z,p]=kcomp(x,f,alpha);
end
if w==1 || w==-1
    J=repmat((1:1:m),m,1);
    I=flipud(rot90(J));
    f=1-abs(I-J)./(m-1); %linear weight
    [po,pe,k,sek,ci,km,ratio,var,z,p]=kcomp(x,f,alpha);
end
if w==2 || w==-1
    J=repmat((1:1:m),m,1);
    I=flipud(rot90(J));
    f=1-((I-J)./(m-1)).^2; %quadratic weight
    [po,pe,k,sek,ci,km,ratio,var,z,p]=kcomp(x,f,alpha);
end



function [po,pe,k,sek,ci,km,ratio,var,z,p]=kcomp(x,f,alpha)
%x: square data matrix
%f: square weight matrix
%alpha: ALPHA - default=0.05.
%po: Observed agreement
%pe: Random agreement
%k: Cohen's kappa
    %k<=0: Poor agreement
    %0.0<k<=0.2: Slight agreement
    %0.2<k<=0.4: Fair agreement
    %0.4<k<=0.6: Moderate agreement
    %0.6<k<=0.8 Substantial agreement
    %0.8<k<=1.0: Perfect agreement
%sekf: kappa error
%ci: kappa confidence interval
%km: Maximum possible kappa, given the observed marginal frequencies
%ratio: kappa observed as proportion of maximum possible
%var: variance
%z: normalized kappa
%p: observed agreement hypotesis test
    %p<0.05: Reject null hypotesis: observed agreement is not accidental
    %p>=0.05: Accept null hypotesis: observed agreement is accidental

m=length(x);
n=sum(x(:)); %Sum of Matrix elements
x=x./n; %proportion
r=sum(x,2); %rows sum
s=sum(x); %columns sum
Ex=r*s; %expected proportion for random agree
pom=sum(min([r';s]));
po=sum(sum(x.*f));
pe=sum(sum(Ex.*f));
k=(po-pe)/(1-pe);
km=(pom-pe)/(1-pe); %maximum possible kappa, given the observed marginal frequencies
ratio=k/km; %observed as proportion of maximum possible
sek=sqrt((po*(1-po))/(n*(1-pe)^2)); %kappa standard error for confidence interval
ci=k+([-1 1].*(abs(-realsqrt(2)*erfcinv(alpha))*sek)); %k confidence interval
wbari=r'*f;
wbarj=s*f;
wbar=repmat(wbari',1,m)+repmat(wbarj,m,1);
a=Ex.*((f-wbar).^2);
var=(sum(a(:))-pe^2)/(n*((1-pe)^2)); %variance
z=k/sqrt(var); %normalized kappa
p=(1-0.5*erfc(-abs(z)/realsqrt(2)))*2;