function polyconvexification
F_ref = [pi/4 0 0 0 0 0 0 0 0]; d = 2; F = F_ref(1:d^2);
r = 4; L = 4; M = 100;
W = @(F)(sum(F.^2,2)-1).^2;
[W_pc,lambda] = multilevel_poly(W,F,r,L,M)

function [W_pc,lambda] = multilevel_poly(W,F,r,L,M)
d = sqrt(size(F,2)); tau_d = (d-1)*d^2+1;
T_F = [F,minors(F)];
delta = r/4; atoms = grid_gen_mat(delta,r,d);
W_pc = 0; lambda = zeros(tau_d,1);
eps_as = 1; ell = 1;
%% Algorithm 9.2
while ell <= L
    W_A = feval(W,atoms);
    T_A = [atoms,minors(atoms)];
    mp = 0;
    while ~mp
        active = active_set(lambda,T_F,eps_as,T_A,W_A,delta,d); % Define N_active
        [lambda,W_pc] = lin_prog(T_F,active,T_A,W_A,tau_d); % Solve optimization problem
        mp = max_princ(lambda,T_F,W_pc,delta,T_A,W_A); % Check maximum principle
        eps_as = eps_as*2; % Increase eps_as if mp not satisfied
    end
    if ell < L
        atoms = refine_coarsen(lambda,T_F,W_pc,T_A,W_A,delta,M,d); % Refinement
        delta = delta/2; eps_as = delta;
    end
    ell = ell+1;
end

function active = active_set(lambda,T_F,eps_as,T_A,W_A,delta,d)
nA = size(T_A,1);
vec = T_A*lambda-W_A;
active = sparse(size(T_A,1),1);
idx_mp = (vec>max(vec)-eps_as);
idx_feas = (max(abs(T_A(:,1:d^2)...
    -ones(nA,1)*T_F(1:d^2)),[],2)<=delta);
active(max(idx_mp,idx_feas)) = 1;

function [lambdas,W_pc] = lin_prog(T_F,active,T_A,W_A,tau_d)
idx = find(active); n_active = nnz(idx);
%% Complete the linear minimization program
...

function mp = max_princ(lambda,T_F,W_pc,delta,T_A,W_A)
vec_A = T_A*lambda-W_A;
mp = ~(max(vec_A)>T_F*lambda-W_pc+delta^2);

function atoms = refine_coarsen(lambda,T_F,W_pc,T_A,W_A,delta,M,d)
vec = T_A*lambda-W_A;
idx = (vec>T_F*lambda-W_pc-M*delta);
atoms = loc_grid_ref_mat(delta/2,T_A(idx,1:d^2),d);

function atoms = grid_gen_mat(delta,r,d)
%% Complete the mesh generation
...

function new_atoms = loc_grid_ref_mat(delta,atoms,d)
nr_old = size(atoms,1); new_atoms = zeros(2^(d^2)*nr_old,d^2);
for z = 0:1:nr_old-1
    ctr = 1;
    for j = 0:1:1
        for k = 0:1:1
            for m = 0:1:1
                for n = 0:1:1
                    new_atoms(z*16+ctr,1) = atoms(z+1,1)+j*delta;
                    new_atoms(z*16+ctr,2) = atoms(z+1,2)+k*delta;
                    new_atoms(z*16+ctr,3) = atoms(z+1,3)+m*delta;
                    new_atoms(z*16+ctr,4) = atoms(z+1,4)+n*delta;
                    ctr = ctr+1;
                end
            end
        end
    end
end

function val = minors(A)
if size(A,2) == 4
    val = A(:,1).*A(:,4)-A(:,2).*A(:,3);
elseif size(A,2) == 9
    val = zeros(size(A,1),10);
    val(:,1:9) = [A(:,5).*A(:,9)-A(:,6).*A(:,8),...
        A(:,4).*A(:,9)-A(:,6).*A(:,7),...
        A(:,4).*A(:,8)-A(:,5).*A(:,7),...
        A(:,2).*A(:,9)-A(:,3).*A(:,8),...
        A(:,1).*A(:,9)-A(:,3).*A(:,7),...
        A(:,1).*A(:,8)-A(:,2).*A(:,7),...
        A(:,2).*A(:,6)-A(:,3).*A(:,5),...
        A(:,1).*A(:,6)-A(:,3).*A(:,4),...
        A(:,1).*A(:,5)-A(:,2).*A(:,4)];
    val(:,10) = A(:,1).*val(:,1)-A(:,2).*val(:,2)...
        +A(:,3).*val(:,3);
end