

function [gamMinCS,gamCNmin,lamdaW_CS,W_CS,W_CN,XminIt,Nit,valFO]= AlgoritmoSDM(x0,alf,A,lamdaA,a)



[n n]=size(A);          % A= matrice Hessiana della Funzione Obiettivo    
W_CS = zeros(2*n,2*n);  % W= matrice del Sistema Lineare Discreto che caratterizza il modello della Rete Neurale associata
                        % all'alg. di BPM 
V=zeros(2*n,1);         % V=V(alf,gam,d) con (W,V) caratterizzano il S. Lineare dinamico
gam1=0;
lamdaW_CS=[];
lamdaW_m=zeros(n,length([0:0.1:1]));
x=0; 
Nit=0;
XminIt = [];
valF0 = 0;

% Condizioni Sufficienti per la convergenza dell'algoritmo: esistenza di gamMinCS che fa si che tutti 
% gli autovalori di W siano complessi, fissato il coefficiente di apprendimento alf:
lamdaA_v = eig(A);        
vettAL = [alf*lamdaA_v]'; 
prodAL = alf*lamdaA;
prodAL_0=[];  
prodAL_1=[]; 
prodAL_2=[];  
gamMvett=[];
     
    for i=1:n,
       gamM=(vettAL(1,i) - 1)^2/(vettAL(1,i) + 1)^2;
       gamMvett=[gamMvett gamM];
    end
    gamMin1=max(gamMvett);
    if (gamMin1 > 0) & (gamMin1 < 1)
        gamMinCS = gamMin1;
    end    
    disp('condizioni sufficienti per la stabilità: il coefficiente momento gam deve essere compreso tra gamMinCS e 1, esclusi, per ogni gam compreso tra ]0,1[');
    %fprintf(1,'\nr', );
    disp('gli autovalori di W corrispondenti sono tutti complessi');
    %fprintf(1,'\nr', ); 
    %calcoliamo adesso gli autovalori di W, per ogni autovalore di A quindi di T
    i=1;
    vettL = [gamMinCS + 0.01: 0.01 : 1]; 
    lamdaW_m=zeros(n,length(vettL));
    for gamMin = gamMinCS + 0.01: 0.01: 1, 
               lamdaT= (1 + gamMin)-(1-gamMin)*alf*lamdaA; 
               y=lamdaT;
               lamdaW1 = (y - sqrt(y^2 - 4*gamMin))/2; 
               lamdaW2 = (y + sqrt(y^2 - 4*gamMin))/2;
               lamdaW = [lamdaW1 lamdaW2];
               lamdaW_m(:,i) = lamdaW.';
               i=i+1;
    end
    lamdaW_CS= lamdaW_m;
    lamdaW_mod=abs(lamdaW_m);
    %in questo caso in cui gli autovalori di W sono tutti complessi, il loro modulo dipende solo dal valore di gamma ammissibile:
    disp('confrontiamo i valori di lamdaW_mod_CS, che è pari alla radice quadrata di gamma');
    disp('con lamdaW_mod calcolato tramite lamdaW, tramite un grafico:');
    lamdaW_mod_CS=sqrt(vettL);
    plot(vettL, lamdaW_mod,'r+',vettL, lamdaW_mod_CS,'b-');
    xlabel('gamma');
    ylabel('modulo lamdaW');
    title('effetti di gamma sul modulo degli autovalori di W:alf=0.1, gamma va scelto tra gamMinCS e 1 ');
    axis([0 1 -1 2]);
    gamMin=gamMinCS + 0.001 ; 
    %calcoliamo la matrice W corripondente a tale gamCNmin
    T= (1 + gamMin)*eye(n,n)-(1-gamMin)*alf*A;
    W1(1:n,:)= [zeros(n,n) eye(n,n)];     %
    W1(n+1:2*n,:)= [-gamMin*eye(n,n)  T];
    W_CS=W1;
    disp('la matrice W corriponde a quella calcolata con il valore del coefficiente gamma uguale a gamMinCS + 0.01');
    %fprintf(1,'\nr', );
    disp('il vettore degli autovalori di W corrisponde ai seguenti valori del coef. gamma: gamMinCS + 0.01: 0.1: 1');
 
    % Consideriamo adesso le Condizioni Necessarie per la convergenza dell'algoritmo che sono 3 
    % la maggiore parte degli autovalori di W sono reali (possono esistere anche autovalori complessi)
                   
         if (prodAL > 0) & (prodAL < 1)     % caso 1 con prodAL = alf*lamdaA 
            x = prodAL;                      
            prodAL_0 = x;
         elseif (prodAL > 1) & (prodAL < 2) % caso 2
            %vettAL(1,i) = x;
            prodAL_1 = prodAL;
         else  %vettAL(1,i)> 2              % caso 3
             prodAL_2 = prodAL; 
             x = prodAL;
             gamCN= (x - 2)/(x + 2);                 % la stabilità del sistema esiste x gamCN < gam < 1
             %vettgam = [vettgam gamCN]';
             %vettgam = gamCN;
         end

         if ~isempty(prodAL_0),
               
                        disp('siamo considerando gli autovalori reali di W (Condizioni Necessarie) dove un autovalore di A é tale che: 0 <alf*lamdaA < 1 (caso 1)');
                        disp('gli autovalori reali di W sono tali che 0 < lamdaW1 < lamdaW2 < 1');
                        disp('il coefficiente momento gam deve essere scelto tra 0 e 1, esclusi');
                        %calcoliamo adesso i due autovalori di W correspondenti a quello di A del caso considerato
                        k1=1;
                        vettL = [0 :0.1: 1]; %!!0.99
                        %lamdaW=zeros(2*n,length(vettL)); 
                        for gamMin = 0: 0.1: 1, %!! 0.99
                            lamdaT= (1 + gamMin)-(1-gamMin)*alf*lamdaA; 
                            y=lamdaT;
                            %ad ogni lambdaT corrispondono 2 lambdaW:
                            lamdaW1 = (y - sqrt(y^2 - 4*gamMin))/2; % con y inteso come y(1,i)
                            lamdaW2 = (y + sqrt(y^2 - 4*gamMin))/2;
                            lamdaW = [lamdaW1 lamdaW2];
                            lamdaW_m(:,k1)= lamdaW.';
                            k1=k1+1;
                        end
                        lamdaW_modCN=abs(lamdaW_m);
                        figure
                        plot(vettL, lamdaW_modCN);
                        xlabel('gamma');
                        ylabel('modulo lamdaW');
                        title('effetti di gamma sul modulo degli autovalori di W: alf = 0.01, gamma va scelto tra 0 e 1, esclusi ');
                        text(vettL(1,length(vettL)-3),lamdaW_m(1,length(vettL)-8), '0 <alf*lamdaA2 < 1');
                        axis([0 1 -1 2]);
                        gamCNmin= 0.01; 
                        %calcoliamo la matrice W corripondente a tale gamCNmin
                        T= (1 + gamCNmin)*eye(n,n)-(1-gamCNmin)*alf*A;
                        W1(1:n,:)= [zeros(n,n) eye(n,n)];     %
                        W1(n+1:2*n,:)= [-gamCNmin*eye(n,n)  T];
                        W_CN=W1;
                        disp('la matrice W corrisponde a quella calcolata con il valore del coefficiente gamma uguale a 0.01');
                        disp('il vettore degli autovalori di W corrisponde ai seguenti valori del coef. gamma: 0: 0.001: 1');
                        % calcoliamo il numero di passi in cui converge l'algoritmo di B.P.M, Nit e il valore della funzione
                        % obiettivo al punto di minimo, valFO
                        %gam = gamCNmin;
                        %[Nit, XminIt, valFO] = SDMom(x0,A,W_CS,alf,a,gam);
                        %fprintf(1,'il nostro algoritmo converge, al punto di minimo XminIt, in un numero di passi pari a Nit = %d \n', Nit); 
        
                        %fprintf(1,'il valore della funzione oggetto, nel punto di minimo, valFO = %d \n', valFO);
                       
          elseif  ~isempty(prodAL_1),   
                        disp('siamo considerando gli autovalori reali di W (C.N.) dove il secondo autovalore di A é tale che: 1 <alf*lamdaA < 2 (caso 2)');
                        %disp('siamo considerando gli autovalori reali di W (C.N.) dove uno degli autovalori di A é tale che: 1 <alf*lamdaA < 2 (caso 2)');
                        disp(' il coefficiente momento gam deve essere scelto tra 0 e 1, esclusi');
                        if prodAl_1 == 1,
                            disp('in questo caso alf*lamdaA è uguale a 1 perciò tutti gli autovalori di W sono complessi');
                        end
                        %calcoliamo adesso gli autovalori di W
                        k2=1;
                        vettL = [0 :0.1: 1]; %!!0.99
                        %lamdaW=zeros(2*n,length(vettL)); 
                        for gamMin = 0: 0.1: 1, %!! 0.99
                            lamdaT= (1 + gamMin)-(1-gamMin)*alf*lamdaA; 
                            y=lamdaT;
                            lamdaW1 = (y - sqrt(y^2 - 4*gamMin))/2; % con y inteso come y(1,i)
                            lamdaW2 = (y + sqrt(y^2 - 4*gamMin))/2;
                            lamdaW = [lamdaW1 lamdaW2];
                            lamdaW_m(:,k2)= lamdaW.';
                            k2=k2+1;
                        end
                        lamdaW_modCN=abs(lamdaW_m);
                        figure
                        plot(vettL, lamdaW_modCN);
                        xlabel('gamma');
                        ylabel('modulo lamdaW');
                        title('effetti di gamma sul modulo degli autovalori di W: alf = 0.2, gamma va scelto tra 0 e 1, esclusi ');
                        text(vettL(1,length(vettL)-3),lamdaW_m(1,length(vettL)-8), '1 <alf*lamdaA2 < 2');
                        axis([0 1 -1 2]);
                        gamCNmin= 0.9; 
                        %calcoliamo la matrice W corripondente a tale gamCNmin
                        T= (1 + gamCNmin)*eye(n,n)-(1-gamCNmin)*alf*A;
                        W1(1:n,:)= [zeros(n,n) eye(n,n)];     %
                        W1(n+1:2*n,:)= [-gamCNmin*eye(n,n)  T];
                        W_CN=W1;
                        disp('la matrice W corriponde a quella calcolata con il valore del coefficiente gamma uguale a 0.9');
                        disp('il vettore degli autovalori di W corrisponde ai seguenti valori del coef. gamma: 0: 0.1: 1');
                        % calcoliamo il numero di passi in cui converge l'algoritmo di B.P.M, Nit e il valore della funzione
                        % obiettivo al punto di minimo, valFO
                        gam = gamCNmin;
                        [Nit, XminIt, valFO] = SDMom(x0,A,W_CN,alf,a,gam);
                       
          else  %prodAL_2 ~= 0 
                         
                       if prodAL_2 ~= 0, %caso 3, ovvero i valori di alf*lamdaA sono mischiati, con prodAL_2 ~= 0)  
                             disp('siamo nel caso delle C.N. con alf*lamdaA > 2 (caso 3);');
                             disp(' il coefficiente momento gam va scelto tra gamCN e 1, esclusi');
                             %gamMin=vettgam(1,1) + 0.01 ;  %tolgliere " +0.01 " ?
                             gamMin=gamCN + 0.01 ;  
                             % NO!!!: T= (1 + gamMin)*eye(n,n)-(1-gamMin)*alf*A;
                             % ma lamdaT direttamente!:
                             lamdaT = (1 + gamMin)-(1-gamMin)*prodAL_2;   % lamdaT associato a prodAL_2(i)        
                             %lamdaT = eig(A);        %vett. colonna
                             %y=lamdaT';  %vett. riga
                             y=lamdaT;
                             ctrl=y^2-4*gamMin;
                             vettL= [gamCN:0.1:1];
                             %vettL= [0:0.1:1];
                             lamdaW_m=zeros(n,length(vettL));
                             if ctrl > 0   % autoval. di W = reali e ad ogni autoval. di A corrisponde 2 autoval. di W
                                %disp('siamo nel caso delle C.N. con alf*lamdaA > 2 (caso 3): il coef. momento gam deve essere scelto tra gamCNmin e 1, esclusi');
                                disp('gli autovalori di W, reali, sono sempre negativi');
                                k1=1;
                                for w=gamCN:0.1:1, 
                                   lamdaW1 = (y - sqrt(y^2 - 4*w))/2;
                                   lamdaW2 = (y + sqrt(y^2 - 4*w))/2;
                                   lamdaW = [lamdaW1 lamdaW2];
                                   lamdaW_m(:,k1)= lamdaW.';
                                   k1=k1+1;
                                end
                                lamdaW_modCN = abs(lamdaW_m);
                                disp('il vettore degli autovalori di W corrisponde ai seguenti valori del coefficiente gamma: gamCNmin: 0.1: 1');
                                %disp('il vettore degli autovalori di W corrisponde ai seguenti valori del coefficiente gamma:0 : 0.1 : 1di cui quelli tra 0 e gamCN sono non ammissibili');
                                figure
                                plot(vettL, lamdaW_modCN);
                                xlabel('gamma');
                                %ylabel('modulo lamdaW');
                                ylabel('lamdaW');
                                %title('effetti di gamma sul modulo degli autovalori di W: alf = 1 perciò gamCNmin < gamma < 1');
                                title('effetti di gamma sugli autovalori di W: alf = 4 perciò gamCNmin < gamma < 1');
                                %title('effetti di gamma sugli autovalori di W:alf = 1 con valori non ammissibili di gam,tra 0 e gamCN');
                                %axis([0 1 -1 3]);
                                axis([0 1 -2 2]);
                                %text(vettL(1,length(vettL)-5),lamdaW_mod(1,length(vettL)-8), 'alf*lamdaA2 > 2');
                                text(vettL(1,length(vettL)-3),lamdaW_m(2,length(vettL)-3), 'alf*lamdaA2 > 2');
                             else    % autovalori di W = complessi
                                disp('siamo nel caso delle C.N. con alf*lamdaA > 2 (caso 3) ma gli autovalori di W sono tutti complessi');
                                disp('in questo caso il modulo degli autovalori di W è uguale alla radice quadrata di gam');
                                %vettL= [0:0.1:1];
                                vettL= [gamCN:0.01:1];
                                k1=1;
                                for w=gamCN:0.01:1, %w=0:0.1:1,
                                   lamdaW1 = (y - sqrt(y^2 - 4*w))/2;
                                   lamdaW2 = (y + sqrt(y^2 - 4*w))/2;
                                   lamdaW = [lamdaW1 lamdaW2];
                                   lamdaW_m(:,k1)= lamdaW.'; %NB:size(lamdaW_m) < a q. dich. all'inizio
                                   k1=k1+1;
                                end
                                lamdaW_modCN = abs(lamdaW_m);
                                disp('il vettore degli autovalori di W corrisponde ai seguenti valori del ...'); 
                                disp('coefficiente gamma: gamCN: 0.01 : 1');
                                figure
                                plot(vettL, lamdaW_modCN);
                                xlabel('gamma');
                                ylabel('modulo lamdaW');
                                %title('effetti coef. momento sul modulo autoval.di W:alf=4 con val. non ammissibili tra 0 e gamCN');
                                %axis([0 1 -2 2]);
                                title('effetti coef. momento sul modulo autoval.di W: alf = 4 perciò gamCNmin=0.84 < gamma < 1');
                               %title('non ci sono autovalori reali ammissibili in questa condizione necessaria');
                                axis([0 1 -1 2]);
                                text(vettL(1,length(vettL)- 8),lamdaW_mod(1,length(vettL)-15), 'alf*lamdaA2 > 2');
                                %NB su 'text': length(vettL) = 16!
                                %text(vettL(1,length(vettL)-4),lamdaW_m(2,length(vettL)-4), 'alf*lamdaA2 > 2');   
                            end % if ctrl
                             %calcoliamo la matrice W corripondente a tale gamCNmin
                             T= (1 + gamMin)*eye(n,n)-(1-gamMin)*alf*A;
                             W1(1:n,:)= [zeros(n,n) eye(n,n)];     %
                             W1(n+1:2*n,:)= [-gamMin*eye(n,n)  T];
                             W_CN=W1;
                             disp('la matrice W_CN corrisponde a quella calcolata con il valore del coefficiente gamma uguale a gamCNmin + 0.01');
                             gamCNmin = gamCN;
                             [Nit,XminIt, valFO] = SDMom(x0,A,W_CN,alf,a,gamMin);
                         end
          end   
       

