function rosenbrockSimplex()
    EPSILON = 0.000001;
    
    [x1, x2] = meshgrid(-2:0.05:2, -1:0.05:2);
    f = rosenbrock(x1, x2);
    figure;
    surfc(x1, x2, f, 'EdgeColor', 'none');
    shading interp;
    view([135, 45]);

    r = x1(ceil(length(x1) * rand));
    sX1 = [r, r, r + 1];
    r = x2(ceil(length(x2) * rand));
    sX2 = [r, r + 1, r];
    sF = [rosenbrock(sX1(1), sX2(1)), rosenbrock(sX1(2), sX2(2)), rosenbrock(sX1(3), sX2(3))];
    S = [sX1' sX2' sF'];
    
    t = 0;
    
    while triArea(S(:, 1), S(:, 2)) > EPSILON
        S = sortrows(S, -3);
        
        sm = sum(S(2:3, 1:2)) / 2;
        
        s1 = 2 * sm - S(1, 1:2);
        y = rosenbrock(s1(1), s1(2));
        
        if y < S(3, 3)
            s1b = 3 * sm - 2 * S(1, 1:2);
            yb = rosenbrock(s1b(1), s1b(2));
            if yb < y
                S(1, :) = [s1b yb];
            else
                S(1, :) = [s1, y];
            end
        else
            s1 = (3 / 2) * sm - 0.5 * S(1, 1:2);
            y = rosenbrock(s1(1), s1(2));
            if y < S(2, 3)
                S(1, :) = [s1, y];
            else
                s1 = 0.5 * (sm + S(1, 1:2));
                y = rosenbrock(s1(1), s1(2));
                if y < S(2, 3)
                    S(1, :) = [s1, y];
                else
                    S(1, 1:2) = 0.5 * (S(1, 1:2) + S(3, 1:2));
                    S(1, 3) = rosenbrock(S(1, 1), S(1, 2));
                    S(2, 1:2) = 0.5 * (S(2, 1:2) + S(3, 1:2));
                    S(2, 3) = rosenbrock(S(2, 1), S(2, 2));
                end
            end
        end
        
        t = t + 1;
        line(S(:, 1), S(:, 2), S(:, 3), 'Color', 'k');
        pause(0.25);
    end
    x = mean(S(:, 1));
    y = mean(S(:, 2));
    z = mean(S(:, 3));
    X = [x, x, x, x, x + 0.025, x - 0.025];
    Y = [y, y, y + 0.025, y - 0.025, y, y];
    Z = [z + 0.025, z - 0.025, z, z, z, z];
    line(X, Y, Z, 'Color', 'y');
    disp(sprintf('Iterations: %d to (%0.5g, %0.5g, %0.5g)\n', t, x, y, z));
    
function v = rosenbrock(x1, x2)
    v = (1 - x1) .^ 2 + 100 * (x2 - x1 .^ 2) .^ 2;
    
function area = triArea(x, y)
    area = abs(x(1) * y(3) - x(1) * y(2) + x(2) * y(1) - x(2) * y(3) + x(3) * y(2) - x(3) * y(1)) / 2;
