https://github.com/sgsellan/developability-of-heightfields
Raw File
Tip revision: 208ddd4aebfd89c943bb863eeb5fef78ab4f4499 authored by Silvia Sellán on 19 June 2021, 18:02:00 UTC
Update README.md
Tip revision: 208ddd4
argmin_X_full.m
function [X,data_2] = argmin_X_full(Z,U,rho,data,A,II,get_energy,usemex,...
    aggregate,opnorm,weights,NN,BXX,BXY,BYX,BYY,BX,BY,C)
% Solves X update in ADMM, which amounts to
% min ||X||_* + (rho/2)||U-AZ-X||^2
% There are many options in this function, but the only one used for
% the results in the paper is aggregate = true, usemex = true, which amounts to just calling
%
%                       AZ = A*Z
%                       X = argmin_X_full_mex(AZ,U,II,rho,weights);
%
%

X = zeros(size(A,1),1);
AZ = A*Z;

if aggregate==1
    if opnorm
        for i=0:(length(II)-1)
            AZi = [AZ(4*i+1),AZ(4*i+2);AZ(4*i+3),AZ(4*i+4)];
            MM = AZi-[U(4*i+1),U(4*i+2);U(4*i+3),U(4*i+4)];
            [T,S,V] = svd(MM);
            h1 = S(1,1);
            if (S(2,2)-(1/rho))>0
                h2 = S(2,2)-(1/rho);
            else
                h2 = 0;
            end
            Hval = T*[h1,0;0,h2]*(V');
            X(4*i+1) = Hval(1,1);
            X(4*i+2) = Hval(1,2);
            X(4*i+3) = Hval(2,1);
            X(4*i+4) = Hval(2,2);
            
        end
    else
        
        if usemex
            
            % weights = ones(length(II),1);
            X = argmin_X_full_mex(AZ,U,II,rho,weights);
            
        else
            
            
            
            for i=0:((size(A,1)/4)-1)
                % disp(i)
                %     yalmip("clear")
                %     H = sdpvar(2,2,'full');
                %     AZi = [AZ(4*i+1);AZ(4*i+2);AZ(4*i+3);AZ(4*i+4)];
                %     objective = norm(H,'nuclear')+(rho/2)*norm([H(1,1);H(1,2);H(2,1);H(2,2)] - ...
                %         AZi+[U(4*i+1);U(4*i+2);U(4*i+3);U(4*i+4)])^2;
                %     options = sdpsettings('verbose',0,'debug',1,'solver','mosek',...
                %         'cachesolvers',1);
                %     constraints = [];
                %     sol = optimize(constraints,objective,options);
                %     Hval = value(H);
                %     X(4*i+1) = Hval(1,1);
                %     X(4*i+2) = Hval(1,2);
                %     X(4*i+3) = Hval(2,1);
                %     X(4*i+4) = Hval(2,2);
                
                AZi = [AZ(4*i+1),AZ(4*i+2);AZ(4*i+3),AZ(4*i+4)];
                MM = AZi-[U(4*i+1),U(4*i+2);U(4*i+3),U(4*i+4)];
                [T,S,V] = svd(MM);
                
                % if S(1,1) == 0
                %     h1 = 0;
                % else
                %     h1 = S(1,1)-(sign(S(1,1))/rho);
                % end
                rho_with_weight = rho/weights(i+1);
                
                
                if S(1,1)-(1/rho_with_weight)>0
                    h1 = S(1,1)-(1/rho_with_weight);
                elseif S(1,1)+(1/rho_with_weight)<0
                    h1 = S(1,1)+(1/rho_with_weight);
                else
                    h1 = 0;
                end
                
                if S(2,2)-(1/rho_with_weight)>0
                    h2 = S(2,2)-(1/rho_with_weight);
                elseif S(2,2)+(1/rho_with_weight)<0
                    h2 = S(2,2)+(1/rho_with_weight);
                else
                    h2 = 0;
                end
                
                Hval = T*[h1,0;0,h2]*(V');
                X(4*i+1) = Hval(1,1);
                X(4*i+2) = Hval(1,2);
                X(4*i+3) = Hval(2,1);
                X(4*i+4) = Hval(2,2);
            end
            
        end
        
    end
elseif aggregate==2
    if opnorm
        for i=0:(length(II)-1)
            AZi = [AZ(4*i+1),AZ(4*i+2);AZ(4*i+3),AZ(4*i+4)];
            MM = AZi-[U(4*i+1),U(4*i+2);U(4*i+3),U(4*i+4)];
            [T,S,V] = svd(MM);
            h1 = S(1,1);
            h2 = (rho/(rho+2))*S(2,2);
            Hval = T*[h1,0;0,h2]*(V');
            X(4*i+1) = Hval(1,1);
            X(4*i+2) = Hval(1,2);
            X(4*i+3) = Hval(2,1);
            X(4*i+4) = Hval(2,2);
            
        end
    else
        
        if usemex
            X = argmin_X_full_L2_mex(AZ,U,II,rho);
        else
            
            for i=0:(length(II)-1)
                
                AZi = [AZ(4*i+1),AZ(4*i+2);AZ(4*i+3),AZ(4*i+4)];
                MM = AZi-[U(4*i+1),U(4*i+2);U(4*i+3),U(4*i+4)];
                [T,S,V] = svd(MM);
                
                if (((rho+2)*S(1,1))-(2*S(2,2)))>0 && (((rho+2)*S(2,2))-(2*S(1,1)))>0
                    h1 = (((rho+2)*S(1,1))-(2*S(2,2)))/(4+rho);
                    h2 = (((rho+2)*S(2,2))-(2*S(1,1)))/(4+rho);
                else
                    assert(S(1,1)>=S(2,2));
                    h1 = (rho/(rho+2))*S(1,1);
                    h2 = 0;
                end
                
                
                
                
                Hval = T*[h1,0;0,h2]*(V');
                X(4*i+1) = Hval(1,1);
                X(4*i+2) = Hval(1,2);
                X(4*i+3) = Hval(2,1);
                X(4*i+4) = Hval(2,2);
                
            end
            
        end
    end
elseif aggregate==3
    for i=0:(length(II)-1)
        
        AZi = [AZ(4*i+1),AZ(4*i+2);AZ(4*i+3),AZ(4*i+4)];
        MM = AZi-[U(4*i+1),U(4*i+2);U(4*i+3),U(4*i+4)];
        [T,S,V] = svd(MM);
        
        
        
        h1 = (rho/(rho+2))*S(1,1);
        h2 = (rho/(rho+2))*S(2,2);
        
        
        
        
        
        Hval = T*[h1,0;0,h2]*(V');
        X(4*i+1) = Hval(1,1);
        X(4*i+2) = Hval(1,2);
        X(4*i+3) = Hval(2,1);
        X(4*i+4) = Hval(2,2);
        
    end
    
end


data_2 = data;
if get_energy
    data_2.energy = 0;
    for i=0:(length(II)-1)
        S = svd([X(4*i+1),X(4*i+2);X(4*i+3),X(4*i+4)]);
        data_2.energy = data_2.energy + sum(S);
    end
end


% for i=1:length(X)
%     % assert(X(i)==Xtest(i))
%     disp(X(i))
%     disp(Xtest(i))
%     pause
% end




end

back to top