function [h,J]=multifilt_gopt(h,Criterion,m_xA,m_xB,r_xxA,r_xxB,w,P_save,options,w_1,o_1,w_2,o_2,w_3,o_3,w_4,o_4,w_5,o_5)
% MULTIFILT_GOPT  Gradient two-group texture optimization
%	
%	[h,J]=multifilt_gopt(h,'Fisher',mx_A,mx_B,r_xxA,r_xxB,w,...
%	    P_save,options) or
%	[h,J]=multifilt_gopt(h,'MeanScatter',mx_A,mx_B,r_xxA,r_xxB,w,...
%	    P_save,options)
%	starts gradient search with the filter h as initial and returns the
%	optimized filter h, and the criterion value J, for either the Fisher
%	or the MeanScatter criterion.. The parameter w is the smoothing 
%	filter to be assumed in the model. m_xA and m_xB are vectors
%	with the mean values of all textures in group A and B, respectively, 
%	and r_xxA and r_xxB are the corresponding autocorrelations stacked
%	on top of each other (all elements assumed to be of equal size).
%	Intermediate results are saved to the file multifilt_gopt_h.mat 
%	with probability P_save at each criterion function evaluation.
%	
%	[h,J]=multifilt_gopt(h,'Criterion',mx_A,mx_B,r_xxA,r_xxB,w,...
%	    P_save,options,w_1,o_1,w_2,o_2,...) 
%	performs the optimization iteratively, first with smoothing filter
%	w_1, then w_2,... and finally with w. Uses options o_1 in first pass,
%	options o_2 in second etc. Useful for speeding up
%	the search, doing fast searches with small smoothing filters
%	initially.
%	
%	For information on the options parameter, see help foptions
%	Additional: options(100)=1 => use constrained search.

% Initialize
[M_h,N_h]=size(h);
[M_w,N_w]=size(w);

% Iterate through list of smoothing filters
i=1;
while exist(sprintf('w_%d',i)) & exist(sprintf('o_%d',i))
  h_=reshape(h.',M_h*N_h,1);
  J=-multifilt_gopt_J(h_,M_h,N_h,Criterion,m_xA,m_xB,r_xxA,r_xxB,w,P_save);
  fprintf(1,'Starting initial gradient search %d: J=%6.4f\n',i,J);
  [h,J]=multifilt_gopt(h,Criterion,m_xA,m_xB,r_xxA,r_xxB,...
    eval(sprintf('w_%d',i)),P_save,eval(sprintf('o_%d',i)));
  i=i+1;
end
h_=reshape(h.',M_h*N_h,1);
J=-multifilt_gopt_J(h_,M_h,N_h,Criterion,m_xA,m_xB,r_xxA,r_xxB,w,P_save);
fprintf(1,'Starting gradient search, J=%6.4f\n',J);


% Parse input arguments
if exist('options')~=1
  options=[];
end
if length(options)>=100,
  constrained=options(100);
  options=trimend(options(1:99));
else
  constrained=0;
end

% Perform the gradient search
if constrained,
  error('Gradient search temporarily disabled due to license shortage');
  % [h_,options]=constr('multifilt_gopt_J',h_,options,[],[],...
  %   'multifilt_gopt_DJ',...
  %   M_h,N_h,Criterion,m_xA,m_xB,r_xxA,r_xxB,w,P_save);
else
  error('Gradient search temporarily disabled due to license shortage');
  % [h_,options]=fminu('multifilt_gopt_J',h_,options,'multifilt_gopt_DJ',...
  %   M_h,N_h,Criterion,m_xA,m_xB,r_xxA,r_xxB,w,P_save);
end

% Compute the final criterion function value
if nargout>1
  J=-multifilt_gopt_J(h_,M_h,N_h,Criterion,m_xA,m_xB,r_xxA,r_xxB,w,P_save);
end

% Transform filter back from vector representation to 2D
h=reshape(h_,M_h,N_h).';


