%  Alzh.m
 %           11 September 2004
 %
 % This is a demo of a simple model of the Alzheimer's disease based on
 % the recurrent Hopfield neural network with autoassociative memory,
 % represented by the weight matrix. The Alzheimer is modelled by "killing
 % synapses",  that is, by setting individual synaptic weights to zero.

 clear
 % definition of triangles
 trnglDefB
 % test plot of the triangle set
 figure(1), clf
 ttl = 'Initial set' ;
 pltSetB(vrt, TB, ttl) 
 
 % matrix TB  (87 x 5) contains specification of 
 %       Ncomp = 5  triangle compositions, i.e, 
 %    five fundamental memories, each described by 87 bits
 m = 87 ;
 Ncomp = 5 ;
 SurvF = 0.9 ;  % survival factor
 
 % encoding fundamental memories
 WW = (TB*TB')/Ncomp ;  %  outer products of fundamental memories
 
 NKilledSynapses = zeros(1, Ncomp); % numbers of killed synapses
 relKill = zeros(1, Ncomp) ;        % ratio of killed synapses
 cntS = zeros(1, Ncomp) ;           % counters will be stored here
 
 for k = 1:Ncomp  % composition loop
    
  % We now gradually destroy the memory matrix by killing a given 
  % proportion of synapses. We multiply each weight by a randomly 
  % generated 0 or 1 

  rKill = round(SurvF*rand(size(WW)));
  W = rKill.*WW;  % a weight matrix with "killed" synapses

  % Now we want to know how many synapses we have killed.
  %   and their fraction
  NKilledSynapses(k) = sum(sum(abs(sign(W-WW)))) ;
  relKill(k) =  NKilledSynapses(k)/(m^2) ;
  
  % The Hopfield relaxation loop should converge to an attractor, 
  % that is, should find a pattern (triangle) stored in a memory 
  % matirx starting from a given initial state.
  % Naturally we may end up in a spurious attractor.
  % The memory matrix is given as W. The initial state is given 
  % by a vector TB(:,k).
  % To know how long the while loop runs we also set a counter.
 
  cnt = 1;
  yTr = sign(W*TB(:,k));  % calculate the output pattern 
                      % for the initial pattern
  dyTr = yTr;         % change in the pattern
  while (max(abs(dyTr))>0) & (cnt < 250)
      yTr22 = sign(W*yTr);
      dyTr = yTr22 - yTr ;
      yTr = yTr22 ;
      cnt = cnt+1 ;
  end

  % the attractors yTr are stored in a matrix Attr
  Y(:, k) = yTr ;

  % the counter is also stored
  cntS(k) = cnt ;

end

% Now we calculate the number of vector elements where 
% the attractors and the compositions differ.
NDiffElements = sum(abs(sign(TB(:,1:Ncomp)- Y))) ;

figure(2), clf
 ttl = 'Final set' ;
 pltSetB(vrt, Y, ttl) 
 
  ['# Killed synapses:  ',  ...
  num2str(NKilledSynapses), sprintf('\n'),  ...
 'Fraction of killed synapses:  ',  ...
  num2str(relKill, '%1.2f '), sprintf('\n'),  ...
 '# Different Elements:  ',  ...
  num2str(NDiffElements), sprintf('\n'),  ...
 '# Relaxation runs:  ',  num2str(cntS)] 
 
% print -f1 -depsc2 AlzhInit
% print -f2 -depsc2 AlzhFinal