function H=lyapunhop(G,C,W,Ib)
%lyapunhop function for Hopfield neural network
% function to plot Lyapunov for Hopfield Neural Net  RWN    02/05/06
% updated 03/06/08 to fix p=2 for n2
% Syntax
%
%   H = lyapunhop(x,y,C,G,W,Ib)
%   info = hoplyapun(code)
%
% Description
%   
%   lyapunov gives a plot of the Lyapunov function for 
%   a 2-vector Hopfield network with tansig(n) activation 
%   functions by calculating
%   -(1/2)a'Wa-Ib'a+sum(Gii*integral(ln(1-x)./(1+x))_betwen_0&a)
%
% typical values from MathCad file: 
% Ib=[0/012446;-0.010841];
% W=[2.15568 -0.0352; -0.0352 2.1522];
% G=2*eye(2); C=5*eye(2);
% n=[0.5;0.25]; a=[0.462117;0.244919];
%
p = input('enter 1 for plot vs amp output a, 2 for plot vs amp input n  ' )
    if isempty(p)
        p=1;
    end
%    
% normalization by C; here C\ gives multiplication by the inverse of C
    Wn=C\W; Gn=C\G; Ibn=C\Ib;
%    
% handle the case of plotting versus activation function output    
 if p==1
      amin=-0.9; ainc=0.1; amax=0.9;
%    amin=-0.8; ainc=0.1; amax=0.8;
    [a1,a2]=meshgrid(amin:ainc:amax);
%
% conversion of a to index: 
% amax=amin+(imax-1)*ainc => imax=1+(amax-amin)/ainc
    imax=1+(amax-amin)/ainc;
% form sum of conductance integrals
    integ1(imax,imax) = Gn(1,1)*quad('log(1-p)./(1+p)',0,a1(imax));
    integ2(imax,imax) = Gn(2,2)*quad('log(1-p)./(1+p)',0,a2(imax));
    integ=integ1+integ2;
%    
% form z = -(0.5)*(a'*W*a)-a'*Ib+integ
z=-(0.5)*(a1.^2*Wn(1,1)+2*a1.*a2*Wn(1,2)+a2.^2*Wn(2,2))-a1*Ibn(1)...
   -a2*Ibn(2)+integ;
%
% form the surface using mesh(x,y,z)
    meshc(a1,a2,z)
    xlabel('a(1)'); ylabel('a(2)'); zlabel('Lyapunov surface');
%    
% handle the case of plotting versus neuron input    
else p==2
    nlim=input('input n1=n2 limits as  [nmin ninc nmax]  ')
    if isempty(nlim)
        nlim=[-5 0.1 5];  
%        nlim=[-2 0.1 2]; 
    end 
        n1=nlim(1):nlim(2):nlim(3);
        n1min=nlim(1); n1inc=nlim(2); n1max=nlim(3);   
%        
        n2=nlim(1):nlim(2):nlim(3);
        n2min=nlim(1); n2inc=nlim(2); n2max=nlim(3); 
% last line above added 03/06/08
    n=[n1 n2]';
    [n1,n2]=meshgrid(n1,n2);
    a1=tanh(n1); a2=tanh(n2);
    imax=1+(n1max-n1min)/n1inc;
    integ1(imax,imax) = Gn(1,1)*quad('log(1-p)./(1+p)',0,a1(imax));
    integ2(imax,imax) = Gn(2,2)*quad('log(1-p)./(1+p)',0,a2(imax));
    integ=integ1+integ2;
%
 z=-(0.5)*(a1.^2*Wn(1,1)+2*a1.*a2*Wn(1,2)+a2.^2*Wn(2,2))-a1*Ibn(1)...
   -a2*Ibn(2)+integ;   
    meshc(n1,n2,z)
    xlabel('n(1)'); ylabel('n(2)'); zlabel('Lyapunov surface');
end
%
% or can also plot usingsurf(x,y,z)
% to view at different angles enter view(az,el) where in degrees
% az is counter clockwise rotation around z axis and el is vertical 
% elevation; view(60,60) gives a good view.