function net = demo_tiny(netCont)
% A demo and test for the sensitivity (d/dx) computation.

% For training we need input vectors and corresponding target vectors.
input  =  [[-1  0  1 1 1 0 -1 -1 0 -.5 .5 .5]; ...
           [-1 -1 -1 0 1 1  1 0 0 -.5 .5 -.5]];
target = [[  0  0  0 0 0 0  0 0 0 -1  -1  1]; ...
          [  1  0 -1 0 0 0  0 0 0  1   1 -1];
          [  0  0  0 0 0 0  0 0 0 -1  -1  1]];

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

% By default, pick random initial weights, and optimize from there:
net = simple_mlp_create([size(input,1) 6 4 size(target,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 = 500;    % Optimization steps between redraws; always restart.
for s = 1:1

  % 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 figure.
  clf();
  visualize2dSurf(net,testx,testy,input,target,1);

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

  pause(0.1);  % 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.




%% Visualize a 2->1 dimensional mapping as a surface.
function visualize2dSurf(net,testx,testy,input,target,ind)
  testInput = [reshape(testx,1,numel(testx)); ...
               reshape(testy,1,numel(testy))];
  testz = simple_mlp_evalMSEr(net, testInput,'feed');
  testz = testz(ind,:);
                                
  % 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;

  figure;
  contour(testx, testy, reshape(testz,size(testx,1),size(testy,2)));
  hold on;

  % Make a quiver plot of gradient field.
  [testx,testy] = meshgrid(-2+.2:.2:2-.2, -2+.2:.2:2-.2);
  testInput = [reshape(testx,1,numel(testx)); ...
               reshape(testy,1,numel(testy))];
  for ii=1:size(testInput,2)
      [dummy,Gx] = simple_mlp_evalMSEr(net, testInput(:,ii), 'feed');
      %disp(Gx)
      testU(ii)=Gx(ind,1);
      testV(ii)=Gx(ind,2);
  end
  quiver (testInput(1,:),testInput(2,:),testU,testV);
%  axis image

%   for ii=1:numel(testx)
%       for jj=1:numel(testx)
%           [dummy,dummy2,u] = simple_mlp_evalMSEr(net, testx(ii),testy(ii));
%       end
%   end
