Matlab solver

input = readmatrix('input_island');

alphas = input(1:1000, 1);
covar = input(1001:end, :);
vals = input(1:1000, 2);
N = length(vals);

betas = optimvar('betas', N, 1);
slack_var = optimvar('slack_var', 1, 'LowerBound', 0);

optimiser = optimproblem('ObjectiveSense', 'maximize');
optimiser.Constraints.maxcons = betas <= alphas;
optimiser.Constraints.mincons = 0 <= betas;
optimiser.Objective = vals' * betas - 100 * slack_var;

optimiser

options = optimoptions(@linprog, 'Display', 'off');
[xLin, fval, exitFlagInt, output] = solve(optimiser,'options',options);

assert(exitFlagInt == 1, "Did not converge");

threshold = 1e-2;
iter = 1;
betas_temp = xLin.betas;
quad_form = betas_temp' * covar * betas_temp;
slack_error = xLin.slack_var;

options = optimoptions(options, 'ConstraintTolerance', 1e-4);
relative_err = 1;

fprintf('Starting loop at %s\n', datestr(now, 'HH:MM:SS'));

while relative_err > threshold
    new_constr = 2 * betas_temp' * covar * betas - slack_var <= betas_temp' * covar * betas_temp;
    new_constr_key = ['iteration', num2str(iter)];
    optimiser.Constraints.(new_constr_key) = new_constr;

    [xLin, fval, exitFlagInt, output] = solve(optimiser, 'options', options);
    assert(exitFlagInt == 1, "Did not converge");

    betas_temp = (betas_temp + xLin.betas) / 2;
    % betas_temp = xLin.betas; % Use this or the one above

    quad_form = xLin.betas' * covar * xLin.betas;
    slack_error = xLin.slack_var;
    iter = iter + 1;

    relative_err = abs((slack_error - quad_form) / quad_form);

    if mod(iter, 100) == 99
        fprintf("Status: %d %d %d %d\n", iter, slack_error, quad_form, relative_err);
    end

    if relative_err <= 0.9
        fprintf("Error reduced <= 90% at iter %d\n", iter);
    end
end

% the above while loop never terminates
fprintf('Ended loop at %s\n', datestr(now, 'HH:MM:SS'));