function simulation3_rebuttal(d, s, p, num_iter)
    spec_gap = 10;
    B = 5;
    normal_sd = 0.1;
    k = 2; min_u1 = 1/(k*sqrt(s));
    rho = 0.1;
    alpha = 5;
    gamma = 0.1;
    tau = 0.1;

    iter = 0;
    exact_recovery_sum = 0; 
    exact_recovery_sum_dt = 0;
    exact_recovery_sum_it = 0;
    exact_recovery_sum_complete = 0;
    while iter < num_iter
        % Construct u1
        % |u_i| = min_u1 + e_i, e_i = a_i/C
        % For sum(|u_i|^2) = 1, C needs to be 
        % ( 2(sum(a_i)) + sqrt(4*(sum(a_i)^2+12*s*sum(a_i^2)) )/( 3*sqrt(s) ).
        a = rand(s, 1);
        C = ( k*sum(a) + sqrt((k^2)*(sum(a)^2)+(k^2)*(k^2-1)*s*sum(a.^2)) )/( (k^2-1)*sqrt(s) );
        e = a/C;
        u1_s = min_u1 + e;
        u1_s = u1_s.*sign(rand(s, 1)-0.5); % Check: sum(u1.^2) == 1?

        % Construct true matrix true_M
        u1 = [u1_s; zeros(d-s, 1)];
        u = [u1, null(u1')*orth(randn(d-1, d-1))];
        L = zeros(d, 1);
        L(2:d) = sort(randn((d-1), 1), 'descend');
        L(1) = L(2) + spec_gap;
        L = diag(L);
        true_M = u*L*u';

        % Construct noise matrix (still without missing) M
        pd = makedist('Normal', 'sigma', normal_sd);
        t = truncate(pd, -B, B);
        noise = triu(random(t,d,d));
        M = true_M + noise + triu(noise,1)';

        % Creating Observation graph
        A = triu(rand(d,d)<=p);
        A = triu(A,1) + A';
        M = M.*A;

        % Optimization
        M_hat = sdp_optim(M, rho, d);
        diag_M_hat = diag(M_hat);
        exact_recovery_sum = exact_recovery_sum + double(sum((diag_M_hat > 1e-4)==[ones(s,1); zeros((d-s),1)])==d);

        [u, v] = TSPCA(M, d, alpha, gamma);

        exact_recovery_sum_dt = exact_recovery_sum_dt + double(sum((abs(u) > 1e-4)==[ones(s,1); zeros((d-s),1)])==d);
        exact_recovery_sum_it = exact_recovery_sum_it + double(sum((abs(v) > 1e-4)==[ones(s,1); zeros((d-s),1)])==d);

        M_complete = matrix_completion_fr(M, A, d, tau);
        M_hat_complete = sdp_optim(M_complete, rho, d);
        diag_M_hat_complete = diag(M_hat_complete);
        exact_recovery_sum_complete = exact_recovery_sum_complete + double(sum((diag_M_hat_complete > 1e-4)==[ones(s,1); zeros((d-s),1)])==d);

        iter = iter+1;
    end
    
    results = [d, s, p, exact_recovery_sum/num_iter, exact_recovery_sum_dt/num_iter, ...
        exact_recovery_sum_it/num_iter, exact_recovery_sum_complete/num_iter];
    
    fileID = fopen(['rebuttal_result.csv'], 'a+');
    fprintf(fileID, '%d, %d, %d, %f, %f, %f, %f\n', results);
    fclose(fileID);
end

function M_hat = sdp_optim(M, rho, d)
    cvx_begin quiet 
        variable M_hat(d,d) hermitian 
        maximize(real(trace(M*M_hat)) - rho*norm(M_hat(:), 1))
        M_hat == hermitian_semidefinite(d);
        trace(M_hat) <= 1;
    cvx_end
end

function [u, v] = TSPCA(M, d, alpha, gamma)
    % Diagonal thresholding SPCA
    alpha_hat = alpha*sqrt(log(d)/d);
    med_diag = median(diag(M'*M)/d);
    J = diag(M'*M)/d >= med_diag*(1+alpha_hat);

    M_reduced = M;
    M_reduced(:,~J) = 0;
    [u,d1] = eigs(M_reduced'*M_reduced, 1);

    rho_hat = max([sum((diag(M'*M)/d).^2-med_diag^2), med_diag^2]);
    tau_hat = sqrt((med_diag*(rho_hat+med_diag))/(d*rho_hat));
    delta_hat = tau_hat*sqrt(2*log(sum(J)));

    u(abs(u)<delta_hat) = 0; 
    
    % Iterative thresholding SPCA
    gamma_hat = gamma*sqrt(max(d1, 1)*log(d)/d);
    diff = 1;
    v0 = u;
    while diff > 1e-4
        [v,~] = svds(M'*M*v0/d, 1);
        v(abs(v)<gamma_hat) = 0;
        diff = norm(v-v0);
        v0 = v;
    end
end

function M_complete = matrix_completion_fr(M, A, d, tau)
    cvx_begin quiet 
        variable M_complete(d,d) hermitian
        minimize(norm((M_complete-M).*A, 'fro')+tau*norm_nuc(M_complete))
    cvx_end
end