function net = demo_spirals(netCont)
% A simple real-time animated demonstration of the two spirals problem.

% For training we need input vectors and corresponding target vectors.
[input,target] = generateTwoSpirals(1,1);
target = .7*target;


% For visualization, we use a more fine-grained grid..
[testx, testy] = meshgrid(-1+.01:.02:1-.01, -1+.01:.02:1-.01);

% By default, pick random initial weights, and optimize from there:
net = simple_mlp_create([2 15 10 1]);
net.regc = 0.011;                    % Regularization coefficient
net.alpha=2;                         % Power of the distance norm
                                     % (for alpha=1 nonsmooth opt needed!)
net.q=2;                             % Distance norm is the Lq norm

if nargin==1, net = netCont; end;    % If there is a parameter,
                                     % retrofit that as a starting point.

uFunc = @simple_mlp_evalMSErVecWrap; % Use the wrapper that hides the 
                                     % structure; weights become a vector.

figure;          % One image that we'll redraw.

atOnce = 200;    % Optimization steps between redraws; always restart.
for s = 1:100

  % Prepare "user parameters" for the optimization algorithms (they are
  % agnostic about the contents of this structure)
  uPar.net = net;
  uPar.inputs = input;
  uPar.targets = target;

  % Around the optimization, we put the weight matrices to a vector, and
  % afterwards we resume the updated weights. This is just a simple
  % concatenation of columns in the matrices; the indices need a bit of 
  % carefulness.
  x = simple_mlp_cell2vec(net.W);
  
  % You can experiment with either the (steepest) Gradient Descent
  % or the Conjugate Gradient method. The latter is *much* faster
  % (but it is to some extent affected by how often it is restarted)
  %x = simple_mlp_optimGD(x, uFunc, uPar, 0.025, atOnce, 0, 0);
  %x = simple_mlp_optimGDh(x, uFunc, uPar, 0.5, atOnce, 0, 0);  
  x = simple_mlp_optimCGP(x, uFunc, uPar, 1, atOnce, 0, 1e-6);
  %x = simple_mlp_optimWrapSolvOpt(x, uFunc, uPar, 1, atOnce, 0, 1e-6);
  
  net.W = simple_mlp_vec2cell(x, net.nneur);

  % The optimization has now progressed for a while. Then we'll just 
  % evaluate some measures to be printed out.
  [f, g] = simple_mlp_evalMSEr(net, input, target);

  netNoReg = net;
  netNoReg.regc = 0;
  fNoReg = simple_mlp_evalMSEr(netNoReg, input, target);
  valREG = f - fNoReg;
  valMSE = fNoReg;
  
  % Update the two-spiral animation.
  clf();
  visualizeSpirals(net,testx,testy,input,target);

  title(sprintf('Iter %d: cost=%f  MSE=%f   Pen.=%f  ',s * atOnce, f, ...
		valMSE, valREG));

  pause(0.2);  % A small pause seems to be needed; otherwise the fig
               % doesn't update... weird but true.. not everything is
               % perfect in Matlab.
end

% ---- End of the main function.



%% This is a demonstration example with some history. 
% You may search the web for "two spirals".
function [input, target] = generateTwoSpirals(density,maxrad)
  indices = 0:96;
  angles = indices * pi / (16*density);
  radii = maxrad * (-indices + (104*density)) / (104*density);
  x = radii .* cos(angles);
  y = radii .* sin(angles);
  input = [[x;y],[-x;-y]];
  target = [-1 * ones(1,97) , 1 * ones(1,97)];


%% Makes a mind-blowing visualization; rather a hack implementation, though.
function visualizeSpirals(net,testx,testy,input,target)
  testz = simple_mlp_evalMSEr(net, [reshape(testx,1,numel(testx)); ...
                                  reshape(testy,1,numel(testy))],'feed');

  % Assume there is a current graphic handle.
  
  % Paint the whole test grid with continuous pixel coloring.
  image(32 + 32 * reshape(testz,size(testx,1),size(testy,2)));
  axis image;
  hold on;

  % Truely "hot" spiral (+1) drawn in red X marks.
  pnow = target > 0;
  plot((0.5+0.5*input(1,pnow))*size(testx,1), ...
       (0.5+0.5*input(2,pnow))*size(testy,2),'rx');
  % "Cold" spiral (-1) drawn in blue marks.
  pnow = target < 0;
  plot((0.5+0.5*input(1,pnow))*size(testx,1), ...
       (0.5+0.5*input(2,pnow))*size(testy,2),'bx');
  
  % Wrong classification -> mark the bad X with a white O.
  trainV =  simple_mlp_evalMSEr(net, input,'feed');
  pnow = (trainV .* target) < 0;
  if sum(pnow)>0
    plot((0.5+0.5*input(1,pnow))*size(testx,1), ...
         (0.5+0.5*input(2,pnow))*size(testy,2),'wo');
  end
