hp_elliptic.m

function hp_model = hp_elliptic( hp_model,model_data )
%function hp_model = hp_elliptic( hp_model,model_data )
%
% compute a h-refinment and rb generation for elliptic problems

      % Generierung des Trainingssets
      rand('state',hp_model.base_model.RB_train_rand_seed);
      %neue Eigenschaft model.training (benötigt für
      %Basisgenerierungsprogramm)
      hp_model.base_model.training = [];
      hp_model.base_model.training = rand_uniform(hp_model.base_model.RB_train_size,...
                           hp_model.base_model.mu_ranges);
                    


      % Erstes eintrag im Baum manuell erstellen
      hp_model.tree = leaf();
      hp_model.tree.anchor = [1;1];   %!!! nicht immer gültig
      hp_model.tree.model =hp_model.base_model;
      hp_model.tree.training = hp_model.base_model.training;
      
            
      littletree = hp_model.tree.anchor;
      
      hp_model.tree = hRB(hp_model.tree);

      % Rekursive funktion zum HP_Baum Aufbau
      % struktur der Anker im Baum werden in littletree (global!)
      % mitgeschrieben um die Funktion part_of_domain zu ermöglichen
      % Diese ist Teil der Traingssetverfeinerung

      function tree_out = hRB(tree_in)
          
      % start with greedy und basisgenerierung
      tree_in.model.training = [tree_in.anchor,tree_in.model.training];
      tree_in.detailed_data = gen_detailed_data(tree_in.model, model_data);


      % Wegen error < errtol
      if  tree_in.detailed_data.RB_info.stopped_on_Nmax == 0
         % error<errtol dann Passt alles 
         
         % Basis und parameter speichern 
         tree_in.parameters = tree_in.detailed_data.RB_info.mu_sequence(:,2:end);
         tree_in.reduced_data = gen_reduced_data(hp_model.base_model,tree_in.detailed_data);
         tree_out = tree_in;
         % Basisgenerierung für Gebiet Abgeschlossen
         disp('--NO  SPLIT--');

      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;


          % Linker Zweig Erbt den Anker vom Elternzweig
          child_left.anchor = tree_in.anchor;
          % Rechte Zweig enthält den ersten neuen greedygefundenen Parameter
            for i=3:child_left.model.RB_train_size
          if ~isequal(child_left.anchor,tree_in.detailed_data.RB_info.mu_sequence(:,i))
             child_right.anchor = tree_in.detailed_data.RB_info.mu_sequence(:,i);
             break;
          end
            end;
              
          % Anker als ersten (wichtig) Eintrag in das neue Trainingsset
          % einfügen
          child_right.model.training =  child_right.anchor;
          child_left.model.training =  child_left.anchor;


            
            
            % 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));
                    distance_r=hp_model.distance_function(child_right.model,child_right.anchor,tree_in.model.training(:,i))  ;
                    % 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.
            littletree=change_littletree(littletree,tree_in.anchor,child_right.anchor);
            
            % Trainingsset verfeinern
            % mittels accept-reject Monte Carlo
          
          %LEFT!
          %boundarybox erstellen:
   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;
          
        
             
          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 ];
              end;
                  
          end;
          
          %RIGHT!
          %boundarybox erstellen:
            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 ];
              end;
                  
          end;
            
            
          
          % Alter Algorithmus zur generierung des Traingssets
          % Sehr schnell aber nicht optimal in der gleichverteilung
          
          % Traingsset verfeinern durch zufälige Punkte im Dreieck aus
          % zufälligen Eckpunkten
          % Nicht Optimal!!!
% 
%                   size_left = length(child_left.model.training(1,:));
%                   size_right = length(child_right.model.training(1,:));
%                   for i=1:(child_left.model.RB_train_size-size_left)
%                         % Zufällige Eckpunkte
%                         rand1 = randi(size_left+i-1,1);
%                         rand2 = randi(size_left+i-1,1);
%                         rand3 = randi(size_left+i-1,1);     
%             
%                         rand4 = rand(1);
%                         rand5 = rand(1)*(1-rand4);
%                         % Neuer Punkt
%                         newpoint = child_left.model.training(:,rand1) + rand4*(child_left.model.training(:,rand2)-child_left.model.training(:,rand1))+rand5*(child_left.model.training(:,rand3)-child_left.model.training(:,rand1));
%                         child_left.model.training = [child_left.model.training,newpoint ];
%                   end;
%             
%                   for i=1:(child_right.model.RB_train_size-size_right)
%                         % Zufällige Eckpunkte
%                         rand1 = randi(size_right+i-1,1);
%                         rand2 = randi(size_right+i-1,1);
%                         rand3 = randi(size_right+i-1,1);     
%             
%                         rand4 = rand(1);
%                         rand5 = rand(1)*(1-rand4);
%                         % Neuer Punkt
%                         newpoint = child_right.model.training(:,rand1) + rand4*(child_right.model.training(:,rand2)-child_right.model.training(:,rand1))+rand5*(child_right.model.training(:,rand3)-child_right.model.training(:,rand1));
%                         child_right.model.training = [child_right.model.training,newpoint ];
%                   end;  

          child_left.training = child_left.model.training;
          child_right.training = child_right.model.training;
        
          %  Rekrusiver Aufruf der neuen Einträge
          %  mit Cells wird ein Baum erzeugt.
          tree_out={hRB(child_left),hRB(child_right)};
      end;
      end


    %Aktuallisiert die Änderungen bei dem Baume mit der Ankerstruktur 
    function tree_out = change_littletree(tree_in,para_expand,para_add)
            t=  iscell(tree_in);
        if t==1
            tree_in{1} = change_littletree(tree_in{1},para_expand,para_add);
             tree_in{2} = change_littletree(tree_in{2},para_expand,para_add);
             tree_out = tree_in;          
        else
            if isequal(tree_in,para_expand)
                tree_out = {tree_in,para_add};
            else
                tree_out = tree_in;
            end;
                
            
        end;
    end

    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 = model.distance_function(model,first_l,param);%norm(first_l-param,2);
            distance_r = model.distance_function(model, first_r,param);%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 - RB Simulation durchführen
           
            if isequal(ltree,anchor)
                is = 1;
            else
                is = 0;
            end;
        end;
        
        
    end

    function value= littletree_first (test)
        t=  iscell(test);
        if t==1
            value=littletree_first(test{1});
        else
            value=test;
        end;
        
    end






end