grid_search.m

function [opt_data,model] = grid_search(model, varargin)
% function opt_data = detailed_grid_search(model, model_data);
%
% Function performing an optimization by grid search using the detailed
% simulations.
%
%Required fields of model:
%model.mu_names: names of the parameters to be optimized
%model.mu_ranges: value ranges for optimization
%model.optimization.init_params: marks which parameters should be optimized
%model.optimization.opt_params.grid_density: Density of the search grid. F.Ex. "3"
%                                means using 3 parameters per dimension
%model.optimitzation.params_to_optimize: Flag vector indicating which
%                           variables to opimize
%
%used variables in Function:
%   P       : cell-array where the used parameters are stored.
%            P{1} = Vector of values for first parameter
%   PM      : Parameter matrix (each row defines a Parameter set for a
%              detailed simulation)
%
%Generated fields of opt_data:
%   output: sequence of outputvalues obtained during optimization
%   parameter_sets: sequence of corresponding parameter sets to output sequence
%   max_output: maximal output value
%   max_output_paramsets: paramset for which the maximal parameter was
%                           obtained
%   
% Markus Dihlmann 03.02.2010



if strcmp(inputname(2),'model_data')
    model_data  = varargin{1};
    detailed_data = [];
        model.optimization.opt_mode = 'detailed';
elseif strcmp(inputname(2),'reduced_data')
    if nargin >1
        reduced_data = varargin{1};
        model_data = varargin{2};
        detailed_data = varargin{3};
            model.optimization.opt_mode = 'reduced';
    else
        reduced_data = varargin{1};
            model.optimization.opt_mode = 'reduced';
    end
end


%Check required fields:
if isempty(model_data)
  model_data = gen_model_data(model);
end;

if isempty(model.mu_names)
    error('No parameter names specified!');
end;

%Total number of parameters defined in model
nbrParams = length(model.mu_names);

%Continue check required fields
if (isempty(model.mu_ranges)||nbrParams~=length(model.mu_ranges))
    error('not the right number of parameter ranges given!');
end;

%if isempty(model.opt_method)
%    error('No optimization method defined!');
%end;

if (isempty(model.optimization.params_to_optimize)||nbrParams~=length(model.optimization.params_to_optimize))
    error('It is not properly defined which parameters are set for optimization and which not!');
end;


% If no initial values are given, set them in the middle of the ranges
if isempty(model.optimization.init_params)
    for i=1:nbrParams
        model.optimization.init_params(i) = (model.mu_ranges{i}(1)+model.mu_ranges{i}(2))/2;
    end;
end;

%Creating search grid
if isempty(model.optimization.opt_params.grid_density)
    model.optimization.opt_params.grid_density =3; %defautl value
end

nbrParams = length(model.mu_names);
nbrOptparams=0; %number of paramters to be optimized

if(size(model.optimization.opt_params.grid_density)==1)
    %same density for all parameter dimensions
    ind=1;
    for i=1:nbrParams
       if model.optimization.params_to_optimize(i)
           stepd = (model.mu_ranges{i}(2)-model.mu_ranges{i}(1))/(model.optimization.opt_params.grid_density-1);
           P{ind} = model.mu_ranges{i}(1):stepd:model.mu_ranges{i}(2);
           ind = ind+1;
       end
    end
    nbrOptParams = ind-1; %number of parameters to be optimized
else
    %for each parameter dimension another density
    ind=1;
    for i=1:nbrParams
        if model.optimization.params_to_optimize(i)
            stepd = (model.mu_ranges{i}(2)-model.mu_ranges{i}(1))/(model.optimization.opt_params.grid_density(i)-1);
           P{ind} = model.mu_ranges{i}(1):stepd:model.mu_ranges{i}(2);
           ind = ind+1; 
        end
    end
    nbrOptParams = ind-1;
end


%Construct ParameterGrid-Matrix
limits=zeros(length(P),1);
for i=1:length(limits)
    limits(i)=length(P{i});
end
nbrCombinations =1;
for i=1:length(limits)
    nbrCombinations=nbrCombinations * limits(i);
end

PM=zeros(nbrCombinations,nbrParams);
counter=ones(1,length(limits)+1);
stop=0;
j=1;
while(~stop)
    ind=1;
    for i=1:nbrParams
        if(model.optimization.params_to_optimize(i))
            PM(j,i)=P{ind}(counter(ind));
            ind = ind+1;
        else
            PM(j,i)=model.optimization.init_params(i);
        end;
    end;
    counter(1)=counter(1)+1;
    for h=1:(length(counter)-1)
              if(counter(h)>limits(h))
                    counter(h)=1;
                    counter(h+1)=counter(h+1)+1;
              end
                
    end
    if counter(length(counter))~=1
                stop=1;
    end
    j=j+1;
end;

%performing a detailed simulation of all parameter sets given by the lines
%of PM
output=zeros(length(PM(:,1)),1);
if strcmp(model.optimization.opt_mode ,'reduced')
    output_der = zeros(length(PM(:,1)),1);
end

for i=1:length(PM(:,1))
    %set_Parameters
    model = model.set_mu(model,PM(i,:));
    
    if(model.verbose>=8)
        disp('Actual parameter set:');
        PM(i,:)
    end;
    
    
    %evaluate output
    if strcmp(model.optimization.opt_mode,'detailed');
        output(i) = model.optimization.objective_function(model,model_data);
    else
        sim_data = lin_evol_rb_derivative_simulation(model, reduced_data);
        sim_data = rb_reconstruction_derivative(model, model_data, detailed_data, sim_data);
        nt=length(sim_data.output);
        output(i) = sim_data.output(nt);
        output_der(i) = sim_data.output_derivative(nt,:) * sim_data.output_derivative(nt,:)';
    end
end;

 
opt_data = [];
opt_data.output = output;
opt_data.parameter_sets = PM;
opt_data.max_output = max(output);
if strcmp(model.optimization.opt_mode ,'reduced')
    opt_data.output_der = output_der;
end

%where to reach this maximum value:
ind=1;
Pmax=[];
ind=search_pos(opt_data.max_output, output,ind)
while(ind>0)
    Pmax=[Pmax;PM(ind,:)];
    ind=search_pos(opt_data.max_output, output,ind+1);
end

opt_data.max_output_paramsets=Pmax;


end

function ind=search_pos(value,vec,start)
    found=0;
    ind=start;
    while ((~found)&&(ind<=length(vec)))
        if(vec(ind)==value)
            found=1;
        else
            ind = ind+1;
        end
    end
    
    if(ind>length(vec))
        ind=-1;
    end
        
end