function simulation2(spec_gap, normal_sd, rho, p, num_iter)
    d = 100;
    s = 50;
    B = 5;
    k = 2; min_u1 = 1/(k*sqrt(s));

    iter = 0;
    exact_recovery_sum = 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);
        
        iter = iter+1;
    end
    results = [spec_gap, normal_sd, rho, p, exact_recovery_sum/num_iter];
    
    fileID = fopen(['simulation2.csv'], 'a+');
    fprintf(fileID, '%f, %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