root/library/utia_legacy/ticket_12/facstr.m @ 1128

Revision 706, 5.6 kB (checked in by smidl, 15 years ago)

eol-native

  • Property svn:eol-style set to native
Line 
1function  [MAPstr, lhs, statistics] = facstr(Fac, Fac0, belief, nbest, ...
2                                             max_nrep, uchns, lambda, order_k)
3% factor-structure estimation
4%
5% [optstr, lhs] = facstr(Fac, Fac0, belief, nbest, max_nrep, uchns, lambda, order_k)
6% [optstr, lhs] = facstr(Fac, Fac0, belief, nbest, max_nrep, uchns, lambda) order_k=2
7% [optstr, lhs] = facstr(Fac, Fac0, belief, nbest, max_nrep) uchns) lambda = 0.9
8% [optstr, lhs] = facstr(Fac, Fac0, belief, nbest, max_nrep) uchns = []
9% [optstr, lhs] = facstr(Fac, Fac0, belief, nbest) max_nrep  = 100 !!!!! change this to 500
10% [optstr, lhs] = facstr(Fac, Fac0, belief) nbest = 1
11% [optstr, lhs] = facstr(Fac, Fac0) belief = 2
12%                               defaults generated by function defaults with the option 's'
13%
14% optstr: maximum a posteriori probability estimate of factor structure
15%%% lhs   : cell vector of the best regressors  -- not true now
16%%%    {1}  values
17%%%    {2}  indicators
18%
19% Fac      : factor          type = 1
20% Fac0     : initial factor  type = 1
21% belief   : user's belief on maximum structure items
22%           (1 items must     be present, 2 items are probably     present
23%            4 items must not be present, 3 items are probably not present)                                 
24% nbest    : how many "best" regressors are maintained
25% max_nrep : maximal number of random starts in search for the best
26%             structure
27% uchns    : list of input channels - if specified, the resulting structure
28%           contains at least one of inputs (suboptimal solution)
29%            Note: channel description can be used instead of "uchns"
30% lambda   : stooping rule threshold
31% order_k  : order of k
32%
33% Design  : L. Tesar. Interface Based on P. Nedoma's previous version of facstrid.
34% Updated : 14.4.2003 - 10.9.2003, MK July 2004
35% Project : post-ProDaCTool
36% Reference: straux1
37   
38% Remark: straux1 uses (or links if it is mex) dydrs and gammaln
39
40%%%% begin patch MK: tests shifted
41[def_belief, def_nbest, def_nrep] = defaults('s');
42
43if nargin<8, order_k  = 2; end;     % LT added stopping rule
44if nargin<7, lambda   = 0.9; end;   % LT added stopping rule
45if nargin<6, uchns  = [];
46elseif iscell(uchns)
47   uu = [];
48   ii = getflds(uchns, 'raction');
49   for i=1:length(ii),
50       if ~isempty(uchns{i}.raction), uu = [uu, uchns{i}.chn]; end
51   end
52   uchns = uu;
53end
54if nargin<5, max_nrep = def_nrep; end;
55if nargin<4, nbest    = def_nbest; end;
56if nargin<3, belief   = def_belief; end;
57if nargin<2, error('facstr needs at least two input parameters'); end;
58
59
60%%%fprintf('FACSTR nruns=%i\n',max_nrep);
61if ~streq(Fac.str,Fac0.str)
62   error('structures of prior and posterior factors in facstrid should equal')
63end
64%%% end patch MK
65
66maxstr = Fac.str;              % structures searched as subselections
67                               % of the current factor structure
68if length(belief)<=1           % belief can be empty
69   if length(belief)==0
70      belief = def_belief;
71   end
72   belief = belief + zeros(1,size(maxstr,2));
73end
74
75npsi  = length(belief);                       % lenght of regressor
76
77%%%% START OF %%%%%%%%%%% LT structure estimation code %%%%%%%%%
78[L, D]   = ld2ld(Fac.LD);
79nu       = Fac.dfm+2; % This is general relation between \nu and dfm
80[L0, D0] = ld2ld(Fac0.LD);
81nu0      = Fac0.dfm+2;
82d        = diag(D);
83d0       = diag(D0);
84
85[optstr1, rgrsout, statistics] = straux1(L, d, nu, L0, d0, nu0, ...
86                                          belief,nbest,max_nrep, lambda, order_k);
87
88
89% result is made as sub-selection of current Factor structure
90maxstr = Fac.str;
91% optstr = maxstr(:,optstr1-1); % probably nobody interested in this now
92
93lhs{1} = [rgrsout.loglik];
94ilh = zeros(length(maxstr),length(rgrsout));
95for f=1:length(rgrsout);
96   ilh(rgrsout(f).strRgr-1,f) = 1;   
97end;
98lhs{2} = ilh;
99%%%% END OF %%%%%%%%%%%%% LT structure estimation code %%%%%%%%%%%
100
101
102% build MAPstr
103   if maxstr(1,npsi)==0,
104      iabs=1;
105      else iabs=0;
106   end
107   ii=find(ilh(1:npsi-iabs,1)==1);
108   if length(ii)~=0,
109       MAPstr=maxstr(:,ii);
110       if (iabs==1 & ilh(npsi,1)==1), MAPstr=[MAPstr [0;1]]; end
111   elseif iabs==1, MAPstr = [0; 1];
112       else MAPstr=[];
113   end
114
115% patch: let the structure contains inputs
116if ~isempty(uchns)
117   [MAPstr, lhs, flag] = checkinp(Fac.str, uchns, lhs);
118end
119
120%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121function [MAPstr, lhs, flag] = checkinp(str, uchns, lhs)
122%  cut best regressors not containing inputs
123%  [MAPstr, lhs, flag] = checkinp(maxstr, uchns, lhs)
124%
125% maxstr: richest factor structure
126% uchns : vector of input channels
127% lhs   : updated information from facstrid
128% flag  : 0 - no inputs
129%         n - 1st ok regressor
130%
131% Design : P. Nedoma
132% Updated: June 2003
133% Project: ProDaCTools remake
134
135vlh   = lhs{1};                        % value of log. likelihood
136ilh   = lhs{2}';                       % the best MAP estimates
137npsi  = size(ilh,1);
138
139% compute MAPstr, set flag
140ii     = find( ilh(1,:) ==1);
141MAPstr = str(:, ii);
142flag   = 0;
143
144nuchn = length(uchns);                 % number of inputs
145iall  = str(1,:);
146ii    = 0*iall;
147
148% build array of indices to input channels in factor structure
149for i=1:nuchn
150    uchn = uchns(i);
151    jj   = find(uchn == iall);
152    if ~isempty(jj), ii(jj) = 0*ii(jj) + 1; end
153end
154
155if all(ii), return; end                % no inputs
156ii = find(ii==1);                      % indices
157
158% find first row containing inputs
159ok = 0;
160for i=1:npsi
161    ind = ilh(i,ii);
162    if any(ind), ok = 1; break; end
163end
164if ~ok, return; end
165
166% cut beginning
167if i>1
168   ilh   = ilh(i:end, :);
169   vlh   = vlh(i:end);
170   flag  = i;
171else
172   flag  = 1; return;
173end
174npsi = size(ilh,1);
175
176% build MAPstr
177ii = find( ilh(1,:) ==1);
178MAPstr = str(:, ii);
179
180% build lhs
181lhs{1} = vlh;
182lhs{2} = ilh';
Note: See TracBrowser for help on using the browser.