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 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)*inc or imax=1+(amax-amin)/inc
imax=1+(amax-amin)/ainc;
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];
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 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.