scm_coercive_greedy.m

Go to the documentation of this file.

function SCM_data=scm_coercive_greedy(model, params, operator,innerproduct_matrix)
%model=scm_greedy(model, model_data)
%
%  SCM greedy function
%  computes the offline data for the calculation of the coercive constant
%  requires model and model data structure
%
%
%   SCM Properties:
%
%
%    SCM_data.Ma                 : number of parameters with an already
%                                  computed constant near to mu
%    SCM_data.Mp                 : number of parameters without computed
%                                  constant near to mu
%    SCM_data.NumberOfParameters : size of the test set used in the SCM
%                                  Greedy procedure
%    SCM_data.eps                : error boundary of the Greedy procedure
%    SCM_data.max_iterations     : maximal number of Greedy iterations
%    SCM_data.parallel           :  = 0  non parallel version
%                                   = 1  parallel (pool stays active after
%                                                               computaton)
%                                   = 2  parallel (closing pool afterwards)
% David Kreplin 05.01.2015

if isempty(params)
    warning('No SCM properties are set -> using standart values')
    %SCM standart values
    params.Ma = inf;
    params.Mp = 0;
    params.NumberOfParameters = 100;
    params.eps=0.1;
    params.max_iterations=inf;
    params.parallel=0;
end


if ~isfield(params, 'Ma')
    params.Ma = inf;
end

if ~isfield(params, 'Mp')
    params.Mp = 0;
end

if ~isfield(params, 'NumberOfParameters')
    warning('SCM: No Number of Parameters is set -> setting NumberOfParameters = 100');
    params.NumberOfParameters = 100;
end

if ~isfield(params, 'eps')
    params.eps = 0.1;
end

if ~isfield(params, 'max_iterations')
    params.max_iterations = inf;
end

if ~isfield(params, 'parallel')
    params.parallel = 0;
end

if ~isfield(params, 'para_rand_seed')
    params.para_rand_seed = 1234;
end

%% Parameter set
if ~isfield(params, 'ParameterSet')
    rand('state',params.para_rand_seed);
    params.ParameterSet = rand_uniform(params.NumberOfParameters,model.mu_ranges);
end
NumberOfParameters = length(params.ParameterSet(1,:));

SCM_data.params = params;

%% Boundary BOX for YLB
model.decomp_mode=1;
[A,~] =  operator(model);
A=remove_dirichlet(A,1,params.dirichlet_gids);

SCM_data.params.Q=length(A);
SCM_data.lower_bound_coercivity = zeros(SCM_data.params.Q,1);
SCM_data.upper_bound_coercivity = zeros(SCM_data.params.Q,1);

%innerproduct_matrix = model.get_inner_product_matrix(model_data);
innerproduct_matrix(SCM_data.params.dirichlet_gids,:)=[];
innerproduct_matrix(:,SCM_data.params.dirichlet_gids)=[];

disp('-- SCM: computing boundary box');
for i = 1:SCM_data.params.Q
    
    [msgstr, msgid] = lastwarn;
    lastwarn('');
    EV= eigs(0.5*(A{1,i}+transpose(A{1,i})),innerproduct_matrix,2,'be');
    [~, warnid] = lastwarn;
    if strcmp(warnid,'MATLAB:eigs:NotAllEigsConverged')
        warning('Eigenvalues did not converged -> using the eig function (more expensive)');
        EV= eig(full(A{1,i}),full(innerproduct_matrix));
    end
    
    lastwarn(msgstr, msgid);
    
    SCM_data.lower_bound_coercivity(i) = min(EV);
    SCM_data.upper_bound_coercivity(i) = max(EV);
end



%% Initalize C
SCM_data.C=[];
SCM_data.C=AddC(params.ParameterSet(:,1),model,operator,innerproduct_matrix,SCM_data);
SCM_data.params.ParameterSet(:,1)=[];
NumberOfParameters = NumberOfParameters-1;
%% SCM Greedy
n=0;
eta=params.eps+1;
%% Greedy Parallel
if params.parallel >= 1
    % s = matlabpool('size');     %OLDER MATLAB VERSION
    % if s == 0
    %     matlabpool('open');
    % end;
    if isempty(gcp('nocreate'))
        parpool;
    end;
    %reduce parallel data amount
    ParameterSet = SCM_data.params.ParameterSet;
    while params.eps <= eta
        n=n+1;
        eta_temp=zeros(NumberOfParameters,1);
        parfor i = 1:NumberOfParameters
            alpha_UB = scm_ub(ParameterSet(:,i),model,operator,SCM_data);
            alpha_LB = scm_coercive_lb(ParameterSet(:,i),model,opereator,SCM_data);
            eta_temp(i) = (alpha_UB-alpha_LB)/(alpha_UB);
        end
        
        [eta,index]=max(eta_temp);
        
        SCM_data.C=AddC(ParameterSet(:,index),model,operator,innerproduct_matrix,SCM_data);
        ParameterSet(:,index) = [];
        NumberOfParameters = NumberOfParameters-1;
        
        disp(['-- SCM: step ' num2str(n) ' - error:' num2str(eta)]);
        if n >= params.max_iterations
            break
        end;
    end
    
    SCM_data.params.ParameterSet=ParameterSet;
    if params.parallel == 2
        %       matlabpool('close');   %OLDER MATLAB VERSION
        delete(gcp('nocreate'))
    end
else
    
    %% non-parallel version
    while params.eps <= eta
        n=n+1;
        eta_temp=zeros(NumberOfParameters,1);
        for i = 1:NumberOfParameters
            alpha_UB = scm_ub(SCM_data.params.ParameterSet(:,i),model,operator,SCM_data);
            alpha_LB = scm_coercive_lb(SCM_data.params.ParameterSet(:,i),model,operator,SCM_data);
            eta_temp(i) = (alpha_UB-alpha_LB)/(alpha_UB);
        end
        
        [eta,index]=max(eta_temp);
        
        SCM_data.C=AddC(SCM_data.params.ParameterSet(:,index),model,operator,innerproduct_matrix,SCM_data);
        
        SCM_data.params.ParameterSet(:,index) = [];
        NumberOfParameters = NumberOfParameters-1;
        
        disp(['-- SCM: step ' num2str(n) ' - error:' num2str(eta)]);
        if n >= SCM_data.params.max_iterations
            break
        end;
    end
    
end



end


function C=AddC(mu,model,operator,innerproduct_matrix,SCM_data)

if isempty(SCM_data.C)
    C.Set = mu;
    model = set_mu(model,mu);
    model.decomp_mode=0;
    [A,~] =  operator(model);
    A=remove_dirichlet(A,0,SCM_data.params.dirichlet_gids);
    A=0.5*(A+transpose(A));
    % [V,D] = eig(full(A),full(innerproduct_matrix));
    try
        [v,d] = fast_eigs( A,innerproduct_matrix,'sm' );
    catch err
        try
            [v,d] = eigs( A,innerproduct_matrix,1,'sm' );
        catch err2
            warning('full spectrum needed, since eigs did not converged');
            [v,d] = eig( full(A),full(innerproduct_matrix) );
            [d,i]=min(max(d));
            v=v(:,i);
        end
    end
    C.coercivity =  d;% min(max(D));
    
    model.decomp_mode=2;
    [A_para,~] = operator(model);
    C.ParameterVector = A_para';
    %v=V(:,end);
    model.decomp_mode=1;
    [A_Matrix,~] =  operator(model);
    A_Matrix=remove_dirichlet(A_Matrix,1,SCM_data.params.dirichlet_gids);
    y_UB=zeros(SCM_data.params.Q,1);
    norm_value = v'*innerproduct_matrix*v;
    for i = 1:SCM_data.params.Q
        y_UB(i) = v'*0.5*(A_Matrix{1,i}+transpose(A_Matrix{1,i}))*v/norm_value;
    end
    C.UBVector = y_UB;
else
    C.Set = [SCM_data.C.Set mu];
    model = set_mu(model,mu);
    model.decomp_mode=0;
    [A,~] =  operator(model);
    A=remove_dirichlet(A,0,SCM_data.params.dirichlet_gids);
     A=0.5*(A+transpose(A));
    %[V,D] = eig(full(A),full(innerproduct_matrix));
    try
        [v,d] = fast_eigs( A,innerproduct_matrix,'sm' );
    catch err
        try
            [v,d] = eigs( A,innerproduct_matrix,1,'sm' );
        catch err2
            warning('full spectrum needed, since eigs did not converged');
            [v,d] = eig( full(A),full(innerproduct_matrix) );
            [d,i]=min(max(d));
            v=v(:,i);
        end
    end
    C.coercivity = [SCM_data.C.coercivity;d];      %min(max(D))];
    
    model.decomp_mode=2;
    [A,~] =  operator(model);
    C.ParameterVector = [SCM_data.C.ParameterVector;A(1:SCM_data.params.Q)'];
    %v=V(:,end);
    model.decomp_mode=1;
    [A_Matrix,~] =  operator(model);
    A_Matrix=remove_dirichlet(A_Matrix,1,SCM_data.params.dirichlet_gids);
    y_UB=zeros(SCM_data.params.Q,1);
    norm_value = v'*innerproduct_matrix*v;
    for i = 1:SCM_data.params.Q
         y_UB(i) = v'*0.5*(A_Matrix{1,i}+transpose(A_Matrix{1,i}))*v/norm_value;
    end
    C.UBVector = [SCM_data.C.UBVector y_UB];
    
end

end

%% Upperbound computation

function alpha_UB = scm_ub(mu,model,operator,SCM_data)

model = set_mu(model,mu);
model.decomp_mode=2;
[theta_q,~] =  operator(model);
alpha_UB=min(theta_q'*SCM_data.C.UBVector);

end


%%

function A=remove_dirichlet(A,mode,I)
if mode == 0
    A(I,:) = [];
    A(:,I)= [];
else
    for i = 1:length(A)
        A{1,i}(I,:) = [];
        A{1,i}(:,I) = [];
    end
end
end