h_refinement_err_ext.m

function hp_model = h_refinement_err_ext( hp_model,model_data )
%function hp_model = h_refinement( hp_model,model_data )
%
% compute the h_refinment of the parameter space 
% extension for error distance
% saves POD base in the little tree  




hp_model.base_model.rho = hp_model.rho;
hp_model.base_model.rho_POD = hp_model.rho_POD;
hp_model.base_model.error_tol1=hp_model.error_tol1;

rand('state',hp_model.base_model.RB_train_rand_seed);
hp_model.base_model.training = rand_uniform(hp_model.base_model.RB_train_size,...
    hp_model.base_model.mu_ranges);

hp_model.tree.anchor = get_mu(hp_model);
hp_model.tree.model = hp_model.base_model;
%hp_model.tree.model.h_part=1;
hp_model.tree.training = hp_model.base_model.training;
%hp_model.tree.model_data = model_data;

%erzeuge basis des anfangsparameters:
training_backup =  hp_model.tree.model.training;
hp_model.tree.model.RB_stop_Nmax = 1;
hp_model.tree.model.training =  hp_model.tree.anchor;

hp_model.tree = create_basis_automatic(hp_model.tree);

hp_model.tree.model.training = training_backup;

% Littletree für die Trainigssetgenerierung
ltreedata.anchor = hp_model.tree.anchor;
ltreedata.detailed_data = hp_model.tree.detailed_data;
ltreedata.reduced_data = hp_model.tree.reduced_data;
ltreedata.model = hp_model.tree.model;
littletree  = ltreedata;



hp_model.tree = hp1(hp_model.tree);






 function tree_out = hp1(tree_in)
        
        %  function tree_out = hp1(tree_in)
        %  generiert die h-Partitionierung für parabollische Probleme rekursiv
        %  tree_in muss schon eine Basis für den Anker besitzen!
        
        %finde den Parameter mit dem größten Fehler -> mu1
        disp('--calculating errors of current domain');
        % Fehler des Trainigssets berechnen:
        post_errs = rb_test_indicator(tree_in.model, ...
            tree_in.detailed_data, tree_in.reduced_data,...
            tree_in.model.training,'');
        
        [max_error, mu_ind]= max(post_errs);
        
        % Fehler von der anchor Approximation (nur relevant bei statischer basisgen):
        error_anchor =  max(rb_test_indicator(tree_in.model, ...
            tree_in.detailed_data, tree_in.reduced_data,...
            tree_in.anchor,...
            ''));
        
        disp(['--Maximum error in current domain: ' num2str(max_error) ' - required: < ' num2str(hp_model.error_tol1)] );
        disp(['--Maximum error of anchor : ' num2str(error_anchor) ' - required: < ' num2str(hp_model.error_tol1/tree_in.model.rho)]);
        
        % Warnen falls Ankerfehler größer als error/rho 
        if error_anchor>hp_model.error_tol1/tree_in.model.rho
            warning('more RB nodes required!');
        end;
        
        if max_error<hp_model.error_tol1
            disp('--NO SPLIT--');
            % fertig, keine zerteilung mehr nötig;
            tree_out = tree_in;
        else
            % Falls error>errortol -> splitten!
            disp('--SPLIT--');
            
            % Erzeuge neue Blätter
            child_right = leaf();
            child_left = leaf();
            % Vererbung von Eigeschaften
            child_right.model = tree_in.model;
            child_left.model = tree_in.model;
         %   child_right.model.h_part =1;
         %   child_left.model.h_part =1;
           % child_right.model_data = tree_in.model_data;
          %  child_left.model_data = tree_in.model_data;
            
            % Linker Zweig erbt den Anker vom Elternzweig
            child_left.anchor = tree_in.anchor;
            % Rechter Zweig enthält den ersten neuen greedygefundenen Parameter
            child_right.anchor = tree_in.model.training(:, mu_ind(1));
            
               % Basis des ersten parameters (mu0) erben
            child_left.detailed_data = tree_in.detailed_data;
            child_left.reduced_data = tree_in.reduced_data;
            
            % Basis für den neuen parameter (mu1) generieren
            backup =  child_right.model.training;
            child_right.model.use_generated_RB_basis = 0;
            
            child_right.model.RB_stop_Nmax = 1;
            child_right.model.training =  child_right.anchor;
            %alternativ: soviel moden nehmen bis error kleiner
            %error_tol/(roh*roh2) mit roh2 >= koerziv(mu)
            
            
            % Basen des Ankers mu1
            child_right = create_basis_automatic(child_right);
            
            child_right.model.training = backup;
            
            child_left.model.training=[];
            child_right.model.training=[];
            
             disp('--dividing trainings set');
            
            % Aufteilen des Elterntrainigsset auf die Kinder
            % hier kommt die d-Funktion zum einsatz            
            leng = length(tree_in.model.training(1,:));
            
            for i = 1:leng
                if ~isequal(tree_in.model.training(:,i),child_right.anchor) && ~isequal(tree_in.model.training(:,i),child_left.anchor)
                    distance_l=hp_model.distance_function(child_left.model,child_left.anchor,tree_in.model.training(:,i),...
                               child_left.detailed_data,child_left.reduced_data);
                    distance_r=hp_model.distance_function(child_right.model,child_right.anchor,tree_in.model.training(:,i),...
                         child_right.detailed_data,child_right.reduced_data)  ;
                    % Aufteilen
                    if distance_l<distance_r
                        child_left.model.training = [child_left.model.training, tree_in.model.training(:,i)];
                    else
                        child_right.model.training = [child_right.model.training, tree_in.model.training(:,i)];
                    end;
                end;
            end;
            
  
            % Protrkolliert die Änderungen in dem kleinen Baum mit.
            ltreedatanew.anchor = child_right.anchor;
            ltreedatanew.detailed_data = child_right.detailed_data;
            ltreedatanew.reduced_data = child_right.reduced_data;
            ltreedatanew.model = child_right.model;
            littletree=change_littletree(littletree,tree_in.anchor,ltreedatanew);
             disp('--refining trainings set (needs some time)');
            
            % Trainingsset mittels accept-reject Monte Carlo verfeinern

            %LEFT SIDE
            % grenzen des Zufallgenerators finden:
            switch hp_model.trainingset_generation_mode
                case 'boundary_box'
                    %boundarybox erstellen:
                    param_min = min(child_left.model.training');
                    param_max = max(child_left.model.training');
                    l = max(param_max-param_min);
                    t = size(child_left.model.training);
                    if isempty(l) || min(l) == 0 || t(2) <= 1
                        box_ranges = child_left.model.mu_ranges;
                    else
                        
                        param_max=param_max+max(l,hp_model.boundarybox_minimum)*hp_model.boundarybox_factor;
                        param_min=param_min-max(l,hp_model.boundarybox_minimum)*hp_model.boundarybox_factor;
                        box_ranges = cell2mat(child_left.model.mu_ranges')';
                        
                        %Box:
                        param_min = max([param_min;box_ranges(1,:)]);
                        param_max = min([param_max;box_ranges(2,:)]);
                        
                        box_ranges =  mat2cell([param_min;param_max]',ones(1,length(param_min)));                       
                    end;
                    
                    
                case 'mu_ranges'
                    box_ranges = child_left.model.mu_ranges;
                otherwise
                    box_ranges = child_left.model.mu_ranges;         
            end;
            
            %Verfeinerung:
            while length(child_left.model.training) < child_left.model.RB_train_size
                
                parameter = rand_uniform(1,box_ranges);
                if part_of_domain(littletree,child_left.anchor,parameter,child_left.model)==1
                    child_left.model.training = [child_left.model.training,parameter ];
                    disp([num2str(length(child_left.model.training)) ' of ' num2str(child_left.model.RB_train_size)]);
                end;
            end;
            
            %RIGHT SIDE
            
            
            switch hp_model.trainingset_generation_mode
                case 'boundary_box'
                    %boundarybox erstellen:
                    param_min = min(child_right.model.training');
                    param_max = max(child_right.model.training');
                    l = max(param_max-param_min);
                    t=size(child_right.model.training);
                    if isempty(l) || min(l) == 0 || t(2) <= 1
                        box_ranges = child_right.model.mu_ranges;
                    else
                        
                        param_max=param_max+max(l,hp_model.boundarybox_minimum)*hp_model.boundarybox_factor;
                        param_min=param_min-max(l,hp_model.boundarybox_minimum)*hp_model.boundarybox_factor;
                        box_ranges = cell2mat(child_right.model.mu_ranges')';
                        
                        %Box:
                        param_min = max([param_min;box_ranges(1,:)]);
                        param_max = min([param_max;box_ranges(2,:)]);
                        
                        box_ranges =  mat2cell([param_min;param_max]',ones(1,length(param_min)));
                        
                    end;
                    
                    
                    
                case 'mu_ranges'
                    box_ranges = child_left.model.mu_ranges;
                otherwise
                    box_ranges = child_left.model.mu_ranges;
                    
            end;
            
            while length(child_right.model.training) < child_right.model.RB_train_size
                
                parameter = rand_uniform(1,box_ranges);
                if part_of_domain(littletree,child_right.anchor,parameter,child_right.model)==1
                    child_right.model.training = [child_right.model.training,parameter ];
                    disp([num2str(length(child_right.model.training)) ' of ' num2str((child_right.model.RB_train_size))]);
                end;
                
            end;
            
         
               
            child_left.training = child_left.model.training;
            child_right.training = child_right.model.training;
            
            
         
            
            % Rekursiver Aufruf
            tree_out={hp1(child_left),hp1(child_right)};
            
        end
        
    end
















    function tree_out = create_basis_automatic(tree_in)
        
        switch (hp_model.find_nr_nodes)
            
            
            case 'automatic'
                
             % generiere basis des ankers -autmatische anzahl der nodes bestimmen;
                % Basisgen. mit POD error:
                RB_ex_algo_backup = tree_in.model.RB_extension_algorithm; %= RB_extension_PCA_fixspace_flexible%
                tree_in.model.RB_extension_algorithm = tree_in.model.RB_extension_algorithm_hp; %=RB_extension_PCA_fixspace_flexible_hp%
                tree_in.detailed_data = gen_detailed_data(tree_in.model,model_data);
                tree_in.reduced_data = gen_reduced_data(tree_in.model,tree_in.detailed_data);
                tree_in.model.RB_extension_algorithm = RB_ex_algo_backup;
                    
                %  Basisgenerierung
                tree_in.model.RB_stop_epsilon = tree_in.model.error_tol1/tree_in.model.rho;
                tree_in.model.RB_stop_Nmax = inf;
                tree_in.model.nr_extension_modes =1;
                % benutze die POD moden:
                tree_in.model.use_generated_RB_basis = 1;

                % Basis genieren (Standart Greedy)
                tree_in.detailed_data = rb_basis_generation(tree_in.model,tree_in.detailed_data);
                %erster Parameter
                tree_in.reduced_data = gen_reduced_data(tree_in.model,tree_in.detailed_data);

                
                %Alte version:
                
                % RB_error:
                
                %    error_RB =  max(rb_test_indicator(tree_in.model, ...
                %        tree_in.detailed_data, tree_in.reduced_data,...
                %        tree_in.anchor,...
                %        ''));

                
                
                %while error_RB > tree_in.model.error_tol1/tree_in.model.rho
                    
                %    tree_in.model.use_generated_RB_basis = 1;
                %    tree_in.detailed_data = rb_basis_generation(tree_in.model,tree_in.detailed_data);
                %    tree_in.reduced_data = gen_reduced_data(tree_in.model,tree_in.detailed_data);
                %    error_RB =  max(rb_test_indicator(tree_in.model, ...
                %        tree_in.detailed_data, tree_in.reduced_data,...
                %        tree_in.anchor,...
                %        ''));
                    
                %end;
                
            case 'static'
                %model.nr_extension_modes werden verwendet - muss gesetzt
                %sein%
                tree_in.detailed_data = gen_detailed_data(tree_in.model,model_data);
                tree_in.reduced_data = gen_reduced_data(tree_in.model,tree_in.detailed_data);
                
                
                
        end;
        
        tree_out = tree_in;
        
    end
    
    


    %Aktuallisiert die Änderungen bei dem Baume mit der Ankerstruktur mit
    %Hilfsfunktionen
    function tree_out = change_littletree(tree_in,para_expand,tree_data_add)
        t=  iscell(tree_in);
        if t==1
            tree_in{1} = change_littletree(tree_in{1},para_expand,tree_data_add);
            tree_in{2} = change_littletree(tree_in{2},para_expand,tree_data_add);
            tree_out = tree_in;
        else
            if isequal(tree_in.anchor,para_expand)
                tree_out = {tree_in,tree_data_add};
            else
                tree_out = tree_in;
            end;
            
            
        end;
    end
    % Überprüft, ob ein Parameter param zu dem Gebiet des anchors gehört 
    % benötigt einen little tree.
    % Hilfsfunktion
    function is = part_of_domain(ltree, anchor,param,model)
        test = iscell(ltree);
        %Navigation
        if test ==1
            % Representanten bestimmen
            first_l=littletree_first(ltree{1});
            first_r=littletree_first(ltree{2});
            % Mit d-Funktion naviagtion durchführen
            distance_l = hp_model.distance_function(model,first_l.anchor,param,...
             first_l.detailed_data,first_l.reduced_data);%norm(first_l-param,2);
            distance_r = hp_model.distance_function(model,first_r.anchor,param,...
                 first_r.detailed_data,first_r.reduced_data);%norm(first_r-param,2);
            % Unterbaum finden und Rekursiv vortfahren
            if distance_l<distance_r
                is = part_of_domain(ltree{1}, anchor,param,model);
            else
                is = part_of_domain(ltree{2}, anchor,param,model);
            end;
        else
            % im passenden Gebiet angekommen
            
            if isequal(ltree.anchor,anchor)
                is = 1;
            else
                is = 0;
            end;
        end;
        
        
    end
    
    % gibt den zugehörigen Anker für die Variable tree zurück.
    function value= littletree_first (tree)
        t=  iscell(tree);
        if t==1
            value=littletree_first(tree{1});
        else
            value=tree;
        end;
        
    end



    
end