opti_lab_3_simp.m

% opti_lab_2.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% tikslo funkcijos ir tt.


% tikslo funkcija
f = @(X) -1 * X(1) * X(2) * X(3);

% lygybinis apribojimas g1 (2xy+2yz+2xz-1)
g1 = @(X) 2*X(1)*X(2) + 2*X(2)*X(3) + 2*X(1)*X(3) - 1;

% kvadratine bauda
b_h1 = @(X) (max(0, -X(1)))^2; % -X(1) yra h1
b_h2 = @(X) (max(0, -X(2)))^2;
b_h3 = @(X) (max(0, -X(3)))^2;
b_g1 = @(X) (g1(X))^2;
kvadratinebauda = @(X) b_h1(X) + b_h2(X) + b_h3(X) + b_g1(X);

% baudos funkcija
baudosf = @(X,r) f(X) + 1/r * kvadratinebauda(X);


% tikslo funkcijos ir tt.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nelder mead simplex

function result = nelderMeadSimplex(point_0, f, r, epsilon)

    result = Result();

    % consts
    alpha = r/100;
    delta_1 = alpha*(sqrt(2+1)+2-1)/(2*sqrt(2));
    delta_2 = alpha*(sqrt(2+1)-1)/(2*sqrt(2));
    coeff_base = 1;
    coeff_expand = 2;
    coeff_shrink = 0.5;
    coeff_shrink_max = -0.5;
    
    % create base simplex
    simplex = [
        point_0(1)          ,     point_0(2)          ,     point_0(3)          ,     NaN; 
        point_0(1) + delta_2,     point_0(2) + delta_1,     point_0(3) + delta_1,     NaN;
        point_0(1) + delta_1,     point_0(2) + delta_2,     point_0(3) + delta_1,     NaN;
        point_0(1) + delta_1,     point_0(2) + delta_1,     point_0(3) + delta_2,     NaN;
    ];
    
    % fill base simplex 4th column with func values
    for i = 1:4
        simplex(i, 4) = f(simplex(i, 1:3), r);
        result.incFuncCallCount();
        result.addIterSimplex(simplex(i, :));
        result.addSimplex(simplex);
    end
    
    % sort descending based on function result
    simplex = sortrows(simplex, -4);

    % function for generating new simplex point
    function point = nextPoint(coeff, simplex)
        centroid = 0.5 * (simplex(3, 1:3) + simplex(4, 1:3));

        % X_naujas = X_h + (1 + coeff)(X_c - X_h)
        new_x_y = simplex(1, 1:3) + (1 + coeff) * (centroid - simplex(1, 1:3));
        point = [new_x_y(1), new_x_y(2), new_x_y(3), f(new_x_y, r)];
        result.incFuncCallCount();
    end

    % iteration loop
    maxiter = 20;
    counter = 0;
    while vecnorm(simplex(4, 1:3) - simplex(3, 1:3)) > epsilon && counter <= maxiter
        counter = counter + 1;
        % create new point
        old_point = simplex(1, :);
        new_point = nextPoint(coeff_base, simplex);

        % adjust new point
        if new_point(4) > simplex(1, 4) % higher than highest
            new_point = nextPoint(coeff_shrink_max, simplex); 

        elseif new_point(4) > simplex(2, 4) % between high and mid
            new_point = nextPoint(coeff_shrink, simplex);

        elseif new_point(4) > simplex(4, 4) % between mid and low
            

        else % lower than lowest
            new_point = nextPoint(coeff_expand, simplex);
        end
        
        simplex(1, :) = new_point;
        result.addIterSimplex(simplex(1, :));
        result.addSimplex(simplex);
        result.addArrow(old_point, new_point);
        simplex = sortrows(simplex, -4);

    end
end

% nelder mead simplex
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% main

start_points = {[0, 0, 0], [1, 1, 1], [0.4, 0.5, 0.2]};
epsilon = 0.0001;
r = 50;

result = nelderMeadSimplex(start_points{1}, baudosf, r, epsilon);
while vecnorm(result.simplexes{end}(4, 1:3) - result.simplexes{end}(3, 1:3)) > epsilon
    r = r/2;
    result = nelderMeadSimplex(result.simplexes{end}(4, 1:3), baudosf, r, epsilon);
end