%MAIN6A        Make a plot of the costfunction surface for 2 weights 

% cvs: $Revision: $

  %%%%%%%%%%%%%%%%%%%%%
  % Initialization    %
  %%%%%%%%%%%%%%%%%%%%%
  
  Ni = 12;                  % Number of external inputs
  Nh = 5;                   % Number of hidden units
  No = 1;                   % Number of output units (classes -1)
  
  % Weight initialization
  range=5;
  
  alpha_i = 1e-5;           % Input weight decay
  alpha_o = 1e-5;           % Output weight decay
  
  iterations = 100; 
  greps = 1e-4;             % Gradient norm stopping criteria

  seed = 45; 
  
  %%%%%%%%%%%%%%%%%%%%
  % Load data        %
  %%%%%%%%%%%%%%%%%%%%
  
  load sp
  x0 = x0(1:12,:)';
  x1 = x1(1:12,:)';
  x2 = x2(1:12,:)';
  y0 = y0';
  y1 = y1';
  y2 = y2';
  

  %%%%%%%%%%%%%%%%%%%%

  [Wi_init,Wo_init] = nr_winit(Ni, Nh, No, range, seed);
  
  [Wi, Wo, Etrain] = nr_trainx(Wi_init, Wo_init, alpha_i, alpha_o, ...
      x0, y0, 0, iterations, greps, 1, 'PR');
  minE = min(Etrain);
  
  %%%%%%%%%%%%%%%%%%%%

  Wi_tmp = Wi;
  Wo_tmp = Wo;
  
  [dummy, index] = sort(-Wi(:));
  w1 = index(2);
  w2 = index(1);
  
  M = 1:30;
  N = 1:30;
  x = logspace(-2,2, length(M)) * Wi(w1);
  y = logspace(-2,1, length(N)) * Wi(w2);
  
  
  cost_c = zeros(length(M),length(N));
  for m = M
    Wi_tmp(w1) = x(m);
    for n = N
      Wi_tmp(w2) = y(n);
      cost_c(m,n) = nr_cost_c(Wi_tmp, Wo_tmp, alpha_i, alpha_o, x0, y0);
    end
    fprintf(1, 'Constructing grid: %d/%d\r', m, max(M))
  end
  
  if 0
    [X, Y] = ndgrid(x,y);
    figure(1)
    contour(X, Y, cost_c)
    xlabel('x')
  end
  
  figure(2)
  %subplot(1,2,1)
  surf(x, y, (cost_c-minE)')
  h = gca;
  set(h, 'XScale', 'log', 'YScale', 'log')
  [i,j] = find(Wi==Wi(w1));
  xlabel(sprintf('Wi(%d,%d)', i, j));
  [i,j] = find(Wi==Wi(w2));
  ylabel(sprintf('Wi(%d,%d)', i, j));
  zlabel('Costfunction - costfunction^{\tau=100}')
  axis([min(x) max(x) min(y) max(y) 0 max(max(cost_c-minE))])
  colormap(cool)
  drawnow
  
  %subplot(1,2,2)
  %contour(x, y, (cost_c-minE)', 50)
  %h = gca;
  %%set(h, 'XScale', 'log', 'YScale', 'log')
  %[i,j] = find(Wi==Wi(w1));
  %xlabel(sprintf('Wi(%d,%d)', i, j));
  %[i,j] = find(Wi==Wi(w2));
  %ylabel(sprintf('Wi(%d,%d)', i, j));
  %zlabel('Costfunction - costfunction^{\tau=100}')
  %colormap(cool)
  %drawnow
  
  
  
  %%%%%%%%%%%%%%%%%%%%%%

    Wi_orig = Wi;
    Wo_orig = Wo;

    i_mask = (Wi == Wi(w1) | Wi == Wi(w2));
    o_mask = (Wo == NaN);

    [Ni,Mi] = size(Wi);
    [No,Mo] = size(Wo);
    D = [Ni,Mi,No,Mo];

    iterations = 11;
    
    w_init = 10^-2;
    
    %%%%%%%%%%%%%%%%%%%%%%

    % Gradient descent

    Wi([w1 w2]) = w_init; 
    clear W_grad1 Cost_grad1
    
    [dWi,dWo] =  nr_gradient(Wi,Wo,alpha_i,alpha_o, x0, y0);
    dWi = dWi .* i_mask;    % Remove non-active terms
    dWo = dWo .* o_mask;

    index = 1;
    while index <= iterations

      W_grad1(index,:) = Wi([w1 w2]);
      Cost_grad1(index) = nr_cost_c(Wi, Wo, alpha_i, alpha_o, x0, y0);
      
      % Determine steplength
      eta = 0.033; 

      % Update weights
      Wi = Wi - eta * dWi;
      Wo = Wo - eta * dWo;

      [dWi,dWo] = nr_gradient(Wi, Wo, alpha_i, alpha_o, x0, y0);
      dWi = dWi .* i_mask;    % Remove non-active terms
      dWo = dWo .* o_mask;
      
      index = index + 1;
    end  


    %%%%%%%%%%%%%%%%%%
    
    Wi = Wi_orig;
    Wo = Wo_orig;

    %%%%%%%%%%%%%%%%%%

    % Plot the "optimization-trace" from the gradient descent
    

    %subplot(1,2,1)
    hold on
    plot3(W_grad1(:,1), W_grad1(:,2), 1.1*(Cost_grad1-minE), 'r--x', ...
	'LineWidth', 3, 'MarkerSize', 15);
    hold off
    %subplot(1,2,2)
    %hold on
    %plot(W_grad1(:,1), W_grad1(:,2), 'r--x', ...
	%'LineWidth', 1, 'MarkerSize', 15);
    %hold off

    %%%%%%%%%%%%%%%%%%%%%%

    % Gradient descent

    Wi([w1 w2]) = w_init; 
    clear W_grad2 Cost_grad2
    
    [dWi,dWo] =  nr_gradient(Wi, Wo, alpha_i, alpha_o, x0, y0);
    dWi = dWi .* i_mask;    % Remove non-active terms
    dWo = dWo .* o_mask;

    index = 1;
    while index <= iterations

      W_grad2(index,:) = Wi([w1 w2]);
      Cost_grad2(index) = nr_cost_c(Wi, Wo, alpha_i, alpha_o, x0, y0);
      
      % Determine steplength
      eta = 0.001; 

      % Update weights
      Wi = Wi - eta * dWi;
      Wo = Wo - eta * dWo;

      [dWi,dWo] = nr_gradient(Wi, Wo, alpha_i, alpha_o, x0, y0);
      dWi = dWi .* i_mask;    % Remove non-active terms
      dWo = dWo .* o_mask;
      
      index = index + 1;
    end  


    %%%%%%%%%%%%%%%%%%
    
    Wi = Wi_orig;
    Wo = Wo_orig;

    %%%%%%%%%%%%%%%%%%

    % Plot the "optimization-trace" from the gradient descent
    

    %subplot(1,2,1)
    hold on
    plot3(W_grad2(:,1), W_grad2(:,2), 1.1*(Cost_grad2-minE), 'g-.d', ...
	'LineWidth', 3, 'MarkerSize', 15);
    hold off
   % subplot(1,2,2)
    %hold on
    %plot(W_grad2(:,1), W_grad2(:,2), 'g-.d', ...
	%'LineWidth', 1, 'MarkerSize', 15);
    %hold off

    %%%%%%%%%%%%%%%%%

    Wi([w1 w2]) = w_init; 
    clear W_pgn Cost_pgn
    
    [dWi,dWo] =  nr_gradient(Wi,Wo,alpha_i,alpha_o, x0, y0);
    dWi = dWi .* i_mask;    % Remove non-active terms
    dWo = dWo .* o_mask;

    
    [dWi,dWo,ddWi,ddWo] =  nr_pseuhess(Wi,Wo,alpha_i,alpha_o,x0,y0);
    dWi = dWi .* i_mask;   % Remove non-active terms
    dWo = dWo .* o_mask;
    
    index = 1;
    while (index <= iterations)

      W_pgn(index,:) = Wi([w1 w2]);
      Cost_pgn(index) = nr_cost_c(Wi, Wo, alpha_i, alpha_o, x0, y0);
      
      % Calculate search direction
      Di = -(dWi ./ ddWi);
      Do = -(dWo ./ ddWo);
      
      % Determine steplength
      eta = nr_linesear(Wi,Wo,Di,Do,alpha_i,alpha_o,x0,y0,12);
      
      % Update weights
      Wi = Wi + eta * Di;
      Wo = Wo + eta * Do;
      
      [dWi,dWo,ddWi,ddWo] = nr_pseuhess(Wi, Wo, alpha_i, alpha_o, x0, y0);
      dWi = dWi .* i_mask;    % Remove non-active terms
      dWo = dWo .* o_mask;
      
      index = index + 1;
    end
    
    %%%%%%%%%%%%%%%%%%
    
    Wi = Wi_orig;
    Wo = Wo_orig;

    %%%%%%%%%%%%%%%%%%

    % Plot the "optimization-trace" from the pseudo_GN
    

    %subplot(1,2,1)
    hold on
    plot3(W_pgn(:,1), W_pgn(:,2), 1.1*(Cost_pgn-minE), 'b-o', 'LineWidth', ...
	3, 'MarkerSize', 10)
    hold off
    %subplot(1,2,2)
    %%hold on
    %plot(W_pgn(:,1), W_pgn(:,2), 'b-o', ...
	%'LineWidth', 1, 'MarkerSize', 15);
    %hold off
    %axis([0 max(x) 0 max(y)])
    
    %%%%%%%%%%%%%%%%%

    set(gcf, 'Renderer', 'zbuffer')
    