vi_detailed_simulation.m

function sim_data = vi_detailed_simulation(model,model_data)
%function sim_data = vi_detailed_simulation(model,model_data)
%
% function performing a detailed simulation of a vi-model (see
% elastic_rope_model) via "quadprog".
%
% Generated fields of sim_data:
%  U: primal solution vector
%  L: Lagrange multiplier vector

% I.Maier 30.05.2011

if ~isfield(model,'use_PDASS')
   model.use_PDASS = 1; 
end

sim_data = [];

model.decomp_mode = 0;
[A,B,f,g] = model.operators(model,model_data);


if ~model.use_PDASS
warning off
options = optimset;
%options.TolFun = 0;
options.MaxIter = 1000; % required for large bases...
options.Display = 'off';
sim_data.U = quadprog(A,-f,B',g,[],[],[],[],[],options);
warning on

sim_data.L = max(0, B\(f - A*sim_data.U) );
else
    [sim_data.U,sim_data.L] = pdass(A, B, f, g);
end





%Primal-Dual-Active-Set_Strategy. Taken from Juliens RBVI Code
function [U,Lambda] = pdass(A, B, f, g)
tol     = 1e-8;
H                       = size(A,1);
index                   = 1:H;
err                     = 1;
U                       = -ones(H,1);
Lambda                  = ones(H,1);
c = 1;
i = 1;

while err>tol
    
    i = i+1;
    Act_set               = index(Lambda+c*(B'*U-g)> 0);    %vorher: ohne c
    Inact_set             = index(Lambda+c*(B'*U-g)<=0);    %vorher: ohne c
    %----Active points----
    U(Act_set)            = -g(Act_set) ;
    Lambda(Inact_set)     = 0 ;
    %---Inactive points---
    M                     = [A(:,Inact_set) , B(:,Act_set)];
    Bg                    = B*g;
    if   size(Act_set,2)>0
      var                 = M\(f-A(:,Act_set)*Bg(Act_set)); ...    % var                 = inv(M)*(f-A(:,Act_set)*Bg(Act_set)); ...
    else        
      var                 = M\f;       %var                 = inv(M)*f; 
    end
    U(Inact_set)          = var(1:length(Inact_set)    );
    Lambda(Act_set)       = var(length(Inact_set)+1:end);
    %---------------------
    err1                  = norm(A*U+B*Lambda-f);
    err2                  = norm(Lambda-max(0,Lambda+c*(B'*U-g)));
    err                   = err1+err2;

    
    if i > 1000
        error('PDASS did not converge within 1000 steps. Please Check');
    end
end