function  [MAPstr, lhs, statistics] = facstr(Fac, Fac0, belief, nbest, ...
                                             max_nrep, uchns, lambda, order_k)
% factor-structure estimation
%
% [optstr, lhs] = facstr(Fac, Fac0, belief, nbest, max_nrep, uchns, lambda, order_k)
% [optstr, lhs] = facstr(Fac, Fac0, belief, nbest, max_nrep, uchns, lambda) order_k=2
% [optstr, lhs] = facstr(Fac, Fac0, belief, nbest, max_nrep) uchns) lambda = 0.9
% [optstr, lhs] = facstr(Fac, Fac0, belief, nbest, max_nrep) uchns = []
% [optstr, lhs] = facstr(Fac, Fac0, belief, nbest) max_nrep  = 100 !!!!! change this to 500
% [optstr, lhs] = facstr(Fac, Fac0, belief) nbest = 1
% [optstr, lhs] = facstr(Fac, Fac0) belief = 2 
%                               defaults generated by function defaults with the option 's' 
%
% optstr: maximum a posteriori probability estimate of factor structure
%%% lhs   : cell vector of the best regressors  -- not true now
%%%    {1}  values
%%%    {2}  indicators
%
% Fac      : factor          type = 1
% Fac0     : initial factor  type = 1
% belief   : user's belief on maximum structure items
%           (1 items must     be present, 2 items are probably     present
%            4 items must not be present, 3 items are probably not present)                                 
% nbest    : how many "best" regressors are maintained
% max_nrep : maximal number of random starts in search for the best
%             structure
% uchns    : list of input channels - if specified, the resulting structure
%           contains at least one of inputs (suboptimal solution)
%            Note: channel description can be used instead of "uchns"
% lambda   : stooping rule threshold
% order_k  : order of k
%
% Design  : L. Tesar. Interface Based on P. Nedoma's previous version of facstrid.
% Updated : 14.4.2003 - 10.9.2003, MK July 2004
% Project : post-ProDaCTool
% Reference: straux1
   
% Remark: straux1 uses (or links if it is mex) dydrs and gammaln

%%%% begin patch MK: tests shifted
[def_belief, def_nbest, def_nrep] = defaults('s');

if nargin<8, order_k  = 2; end;     % LT added stopping rule
if nargin<7, lambda   = 0.9; end;   % LT added stopping rule
if nargin<6, uchns  = [];
elseif iscell(uchns)
   uu = [];
   ii = getflds(uchns, 'raction');
   for i=1:length(ii),
       if ~isempty(uchns{i}.raction), uu = [uu, uchns{i}.chn]; end
   end
   uchns = uu;
end
if nargin<5, max_nrep = def_nrep; end;
if nargin<4, nbest    = def_nbest; end;
if nargin<3, belief   = def_belief; end;
if nargin<2, error('facstr needs at least two input parameters'); end;


%%%fprintf('FACSTR nruns=%i\n',max_nrep);
if ~streq(Fac.str,Fac0.str)
   error('structures of prior and posterior factors in facstrid should equal')
end
%%% end patch MK

maxstr = Fac.str;              % structures searched as subselections
                               % of the current factor structure
if length(belief)<=1           % belief can be empty
   if length(belief)==0
      belief = def_belief; 
   end
   belief = belief + zeros(1,size(maxstr,2));
end

npsi  = length(belief);                       % lenght of regressor

%%%% START OF %%%%%%%%%%% LT structure estimation code %%%%%%%%%
[L, D]   = ld2ld(Fac.LD);
nu       = Fac.dfm+2; % This is general relation between \nu and dfm
[L0, D0] = ld2ld(Fac0.LD);
nu0      = Fac0.dfm+2;
d        = diag(D);
d0       = diag(D0);

[optstr1, rgrsout, statistics] = straux1(L, d, nu, L0, d0, nu0, ...
                                          belief,nbest,max_nrep, lambda, order_k);


% result is made as sub-selection of current Factor structure
maxstr = Fac.str;
% optstr = maxstr(:,optstr1-1); % probably nobody interested in this now

lhs{1} = [rgrsout.loglik];
ilh = zeros(length(maxstr),length(rgrsout));
for f=1:length(rgrsout);
   ilh(rgrsout(f).strRgr-1,f) = 1;   
end;
lhs{2} = ilh;
%%%% END OF %%%%%%%%%%%%% LT structure estimation code %%%%%%%%%%%


% build MAPstr
   if maxstr(1,npsi)==0,
      iabs=1;
      else iabs=0;
   end
   ii=find(ilh(1:npsi-iabs,1)==1);
   if length(ii)~=0,
       MAPstr=maxstr(:,ii);
       if (iabs==1 & ilh(npsi,1)==1), MAPstr=[MAPstr [0;1]]; end
   elseif iabs==1, MAPstr = [0; 1];
       else MAPstr=[];
   end

% patch: let the structure contains inputs
if ~isempty(uchns)
   [MAPstr, lhs, flag] = checkinp(Fac.str, uchns, lhs);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [MAPstr, lhs, flag] = checkinp(str, uchns, lhs)
%  cut best regressors not containing inputs
%  [MAPstr, lhs, flag] = checkinp(maxstr, uchns, lhs)
%
% maxstr: richest factor structure
% uchns : vector of input channels
% lhs   : updated information from facstrid
% flag  : 0 - no inputs
%         n - 1st ok regressor
%
% Design : P. Nedoma
% Updated: June 2003
% Project: ProDaCTools remake

vlh   = lhs{1};                        % value of log. likelihood
ilh   = lhs{2}';                       % the best MAP estimates
npsi  = size(ilh,1);

% compute MAPstr, set flag
ii     = find( ilh(1,:) ==1);
MAPstr = str(:, ii);
flag   = 0;

nuchn = length(uchns);                 % number of inputs
iall  = str(1,:);
ii    = 0*iall;

% build array of indices to input channels in factor structure
for i=1:nuchn
    uchn = uchns(i);
    jj   = find(uchn == iall);
    if ~isempty(jj), ii(jj) = 0*ii(jj) + 1; end
end

if all(ii), return; end                % no inputs
ii = find(ii==1);                      % indices

% find first row containing inputs
ok = 0;
for i=1:npsi
    ind = ilh(i,ii);
    if any(ind), ok = 1; break; end
end
if ~ok, return; end

% cut beginning
if i>1
   ilh   = ilh(i:end, :);
   vlh   = vlh(i:end);
   flag  = i;
else
   flag  = 1; return;
end
npsi = size(ilh,1);

% build MAPstr
ii = find( ilh(1,:) ==1);
MAPstr = str(:, ii);

% build lhs
lhs{1} = vlh;
lhs{2} = ilh';