function niemdiss_pre(setup)
% This function (re-)creates the simulated datasets used in the thesis
% illustrations, and writes them in the datadir. The random seed is
% fixed for each case so this should be exactly repeatable at least
% with the same software and hardware installation.

%% Generate a linearly separable two-dimensional case:
rand('seed',123);
par.indim = 2;
iswp.x = [.25,.25 ];
iswp.n = [2,1     ];
iswp.sep = 0.1; % separated by at least 0.2 at each point.
d = generate_datasets([200,200], @gener_box2, par, @isw_2dvshape, iswp);

figure;
vis_binclass_dataset(d);
axis([-1 1 -1 1],'square')

dlmwrite(fullfile(setup.datadir,'2dlinsep.data'),d,' ');

%% Non-linear, but separable two-dimensional case ("sharp V-shape"):
rand('seed',123);
par.indim=2;
iswp.x = [0,-.5 ; 0,-.5 ];
iswp.n = [1,1   ; -1,1  ];
iswp.sep = 0.1; % separated by at least 0.2 at each point.
d = generate_datasets([200,200], @gener_box2, par, @isw_2dvshape, iswp);

figure;
vis_binclass_dataset(d);
axis([-1 1 -1 1],'square')

dlmwrite(fullfile(setup.datadir,'2dv.data'),d,' ');

%% Non-linear, separable, 2D, requires many neurons. ("sharp W-shape"):

rand('seed', 116)
par.indim=2;
iswp.x1 = [-.5,-.5 ; -.5,-.5 ];
iswp.x2 = [.5,-.5 ; .5,-.5 ];
iswp.n1 = [2,1   ; -2,1    ];
iswp.n2 = [2,1   ; -2,1     ];
iswp.sep = 0.1; % separated by at least 0.2 at each point.
d = generate_datasets([200,200], @gener_box2, par, @isw_2dwshape, iswp);

figure;
vis_binclass_dataset(d);
axis([-1 1 -1 1],'square')

dlmwrite(fullfile(setup.datadir,'2dw.data'),d,' ');

%% Generate a nonseparable, sharp V-shape in two dimensions:
% This is used in the thesis Intro to illustrate nonseparable
% classes (2 neurons) and overfitting (10 neurons)
rand('seed', 125)
par.indim=2;

iswp.x = [0,-.5 ; 0,-.5 ];
iswp.n = [1,1   ; -1,1     ];
iswp.sep = -0.2; % nonseparable.
d = generate_datasets([100,100], @gener_box2, par, @isw_2dvshape, iswp);

figure;
vis_binclass_dataset(d);
axis([-1 1 -1 1],'square')

dlmwrite(fullfile(setup.datadir,'2dvOverlap.data'),d,' ');

%% 2D, 3-class, no overlap, but imbalance and class noise in all classes
rand('seed',128);
par.indim=2;
d = generate_datasets([20,80,300], @gener_box2, par, @isw_easyballs);
d(:,3) = classnoise(d(:,3),[.2,.2,.2]);

figure;
vis_binclass_dataset(d);
axis([-1 1 -1 1],'square')

dlmwrite(fullfile(setup.datadir,'2d3cClassnoise.data'),d,' ');


%% 2D, 3-class, overlap and imbalance
rand('seed',128);
par.indim=2;
d = generate_datasets([10,40,950], @gener_box2, par, @isw_balls);

figure;
vis_binclass_dataset(d);
axis([-1 1 -1 1],'square')

dlmwrite(fullfile(setup.datadir,'2d3cOverlapImbalanced.data'),d,' ');

%% 2D, 3-class, overlap, imbalance, and class noise
rand('seed',128);
par.indim=2;
d = generate_datasets([10,40,950], @gener_box2, par, @isw_balls);
d(:,3) = classnoise(d(:,3),[.0,.05,.05]);

figure;
vis_binclass_dataset(d);
axis([-1 1 -1 1],'square')

dlmwrite(fullfile(setup.datadir,'2d3cOverlapImbalancedClassnoise.data'),d,' ');

%% Similar, but not identical, to the earlier case study (ESANN2016):
rand('seed',128);
par.indim=2;
d = generate_datasets([20,40,340], @gener_box2, par, @isw_esann16);

figure;
vis_binclass_dataset(d);
axis([-1 1 -1 1],'square')

dlmwrite(fullfile(setup.datadir,'2d3cOverlapImbalanced400.data'),d,' ');

%% Small data for initial, fast computations.:
rand('seed',128);
par.indim=2;
d = generate_datasets([5,10,85], @gener_box2, par, @isw_esann16);

figure;
vis_binclass_dataset(d);
axis([-1 1 -1 1],'square')

dlmwrite(fullfile(setup.datadir,'2d3cOverlapImbalanced100.data'),d,' ');



%% Helper functions ------------

function b = isw_2dwshape(p, ic, par)
  % Combine two "V"-shapes to make a "W".
  % (Works for two convex polyhedra)
  %
  % par.x1 = one point belonging to each hyperplane as row mat
  % par.x2 = one point belonging to each hyperplane as row mat
  %         
  % par.n1 = directions of the hyperplanes as row matrix
  % par.n2 = directions of the hyperplanes as row matrix
  %
  % par.sep = separability (>0 means separable; <0 overlap)

  par1.x = par.x1;
  par1.n = par.n1;
  par1.sep = par.sep;
  par2.x = par.x2;  
  par2.n = par.n2;
  par2.sep = par.sep;

  if (ic==1) % within either of two shapes:
    b = isw_2dvshape(p, ic, par1) || isw_2dvshape(p, ic, par2);
  elseif (ic==2) % within both "outside shapes"..
    b = isw_2dvshape(p, ic, par1) && isw_2dvshape(p, ic, par2);
  else
    error('This works only for 2 classes. Requested class %d', ic);
  end




function b = isw_2dvshape(p, ic, par)
  % Essentially works for any convex polyhedron defined as
  % hyperplanes
  %
  % par.x = one point belonging to each hyperplane as row mat
  %         
  % par.n = directions of the hyperplanes as row matrix
  %         (will be re-normalized here for convenience)
  %
  % par.sep = separability (>0 means separable; <0 overlap)

  indim = size(par.x,2);
  nplanes = size(par.x,1);
  assert(size(par.n,1) == nplanes, 'dimension mismatch');
  assert(size(par.n,2) == indim, 'dimension mismatch');
  
  sep = par.sep;
  % FIXME: Doubt this implementation before verification
  % (seems to work as of now, though)
  
  % re-normalize:
  n = par.n ./ repmat(sqrt(sum(par.n .* par.n,2)),1,indim);
  % translate points:
  tp = repmat(p,nplanes,1) - par.x;
  % compute dot prod of the translated point with each direction:
  dotps = sum(n .* tp, 2);

  if (ic==1) % one side of all planes, thereby "within polyhedron"
    b = sum(dotps >= (0+sep)) == nplanes;
  elseif (ic==2) % other side of any plane, thereby "outside"
    b = sum(dotps < (0-sep)) > 0;
  else
    error('This works only for 2 classes. Requested class %d', ic);
  end
  
% Old version with fixed two planes:
%  n1 = par.n1 / norm(par.n1,2);
%  n2 = par.n2 / norm(par.n2,2);
%  sep = par.sep;
  
%  dist1 = (p - par.x1) * (n1');
%  dist2 = (p - par.x2) * (n2');

function iswithin = isw_esann16(p,ic,par)
  % Parameters unused here; could have the ball parameters there.
  if ic==1
    iswithin = within_ball(p, [0.25 0.25], 0.5);
  elseif ic==2
    iswithin = within_ball(p, [-0.25 -0.25], 0.5);
  elseif ic==3
    iswithin = 1;
  else
    error('This function is for classes 1,2,3. Class %d requested.', ic);
  end

  
function iswithin = isw_easyballs(p,ic,par)
  % Parameters unused here; could have the ball parameters there.
  if ic==1
    iswithin = within_ball(p, [-.5 0.0], 0.3);
  elseif ic==2
    iswithin = within_ball(p, [0.25 0.5], 0.4);
  elseif ic==3
    iswithin = within_ball(p, [0.5 -0.5], 0.5);
  else
    error('This function is for classes 1,2,3. Class %d requested.', ic);
  end


function iswithin = isw_balls(p,ic,par)
  % Parameters unused here; could have the ball parameters there.
  if ic==1
    iswithin = within_ball(p, [0.1 0.1], 0.1);
  elseif ic==2
    iswithin = within_ball(p, [-0.25 -0.25], 0.5) ...
	|  within_ball(p, [0.5 -0.3], 0.2) ...
	|  within_ball(p, [0.5 0.5], 0.2);
  elseif ic==3
    iswithin = 1;
    iswithin = ~within_ball(p, [0.1 0.1], 0.1) ...
	& ~within_ball(p, [-0.25 -0.25], 0.5);
  else
    error('This function is for classes 1,2,3. Class %d requested.', ic);
  end
  
function binsel = within_ball(d, c, r)
  x = d - repmat(c, size(d,1), 1);
  binsel = sum(x .* x, 2) < r^2;
return




% Unused, as of now.. we do all computations in the main C++ implementation.
function mlp = UNUSEDtrain_simple(nneur, input, target)
  iters =  10000; %max iterations
  atOnce = 42; %line searches before restart
  net = simple_mlp_create(nneur);
  net.regc = 0.00005;                  % 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

  for ilay=1:(numel(nneur)-2)
    net.activ{ilay}  = @simple_mlp_logsigfwd; % TODO: Check logsig
                                              % codes once more
    net.activD{ilay} = @simple_mlp_logsigderiv;
  end

  % (Output layer already has linear activation by default.)

  uFunc = @simple_mlp_evalMSErVecWrap; % Wrapper hiding structure
  uPar.net = net;
  uPar.inputs = input;
  uPar.targets = target;

  besterr = Inf;
  mlp = net;
  for ii=1:1  % multistart it
    disp(ii);
    x = simple_mlp_cell2vec(net.W);
    x = -.25 + .5*rand(size(x));
    for iit = 1:(iters/atOnce)
      x = simple_mlp_optimCGP(x, uFunc, uPar, 1, atOnce, 0, 1e-9);
    end

    net.W = simple_mlp_vec2cell(x, net.nneur);
    err = simple_mlp_evalMSEr(net, input, target);
    if (err < besterr) mlp = net; end
  end

