%% RBFN update: new points adding

ga=0;
gr=0;
rnd=0;

%% adding solution of aprox optimization point
points(npoints+1,:) = ga_solution(:);
xreal(:) = points(npoints+1,:)'.*(Domain(:,2)-Domain(:,1))+Domain(:,1);
ptsvals(npoints+1) = F(xreal(:));
if ptsvals(npoints+1)<realbsf_value
    realbsf(:) = points(npoints+1,:);
    realbsf_value = ptsvals(npoints+1);
%     Feval(npoints, realbsf_value, realbsf);
end
npoints = npoints+1;
ga = 1;
ga1 = ga2;
ga2 = npoints;

%% adding point in the direction of the better of two last optim. solutions
if ga
    dif = points(ga1) - points(ga2);
    if ptsvals(ga1) < ptsvals(ga2)
        points(npoints+1,:) = points(ga1,:) + dif(:);
    else
        points(npoints+1,:) = points(ga2,:) - dif(:);
    end
    if Return_to_domain
        for i=1:Dim
            if points(npoints+1,i) > 1 points(npoints+1,i) = 1; end
            if points(npoints+1,i) < 0 points(npoints+1,i) = 0; end
        end
    end
        
    xreal(:) = points(npoints+1,:)'.*(Domain(:,2)-Domain(:,1))+Domain(:,1);
    ptsvals(npoints+1) = F(xreal(:));
    if ptsvals(npoints+1)<realbsf_value
        realbsf(:) = points(npoints+1,:);
        realbsf_value = ptsvals(npoints+1);
%         Feval(npoints, realbsf_value, realbsf);
    end
    npoints = npoints+1;
    gr = 1;
    grid = npoints;
end

%% adding random point
points(npoints+1,1:Dim) = rand(1,Dim);
xreal(:) = points(npoints+1,:)'.*(Domain(:,2)-Domain(:,1))+Domain(:,1);
ptsvals(npoints+1) = F(xreal(:));
if ptsvals(npoints+1)<realbsf_value
    realbsf(:) = points(npoints+1,:);
    realbsf_value = ptsvals(npoints+1);
%     Feval(npoints, realbsf_value, realbsf);
end
npoints = npoints+1;
rnd = 1;
rndid = npoints;

%% display results
if rbfndisplayconf
    news = sprintf( '%12d   %10g', npoints, realbsf_value);
    if ga 
        news = sprintf( '%s   %10g', news, ptsvals(ga2) );
    else
        news = sprintf( '%s             ', news);
    end
    if gr 
        news = sprintf( '%s   %10g', news, ptsvals(grid) );
    else
        news = sprintf( '%s             ', news);
    end
    if rnd 
        news = sprintf( '%s   %10g', news, ptsvals(rndid) );
    else
        news = sprintf( '%s             ', news);
    end
    news = sprintf( '%s   %10g', news, rbfrad );
    display( sprintf( '%s', news) );
end
    
%% RBFN plot

if showrbfn==1
    switch Dim
        case 1
            N=100;
            aval = zeros(N+1,1);
            x = 0:1/N:1;
            for i = 1:N+1
                aval(i) = Avalue(x(i));
            end
            plot(x(:),aval(:),points(1:npoints-3,1),ptsvals(1:npoints-3),'+r',points(npoints-2:npoints,1),ptsvals(npoints-2:npoints),'*k');
            if showrbfnsave==1
                filename = sprintf( 'rbfn_%02d.png', istep );
                saveas( gcf, filename, 'png' );
                filename = sprintf( 'rbfn_%02d.eps', istep );
                saveas( gcf, filename, 'eps' );
            end
            pause
        case 2
            N=100;
            aval = zeros(N+1);
            x = 0:1/N:1;
            [X,Y] = meshgrid(x);
            for i = 1:N+1
                for j = 1:N+1
                    aval(i,j) = Avalue([(i-1)/N (j-1)/N]);
                end
            end
            mesh(X,Y,aval);
            pause
        otherwise
            display( 'Dimension of problem is too high for ploting of approximation' );
    end
end

%% author: Anicka Kucerova, Czech Technical University in Prague
%% version 0.1.2
%% 30 Mar 2009