001 function [X,res,comb,f]=rfc2load_fat(f,res,num_cc)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034 ni = nargin;
035 no = nargout;
036 error(nargchk(1,3,ni));
037
038 if ni<2, res=[]; end
039 if ni<3, num_cc=[]; end
040
041
042 if ~(round(f)==f)
043 warning('The frequency matrix is not interger matrix. The matrix will be scaled.')
044 f=floor(100*f);
045 end
046
047 N=length(f);
048
049
050
051 if isempty(res)
052 res0=fr2res(f);
053
054
055 [RFC,RFC1,res] = tp2rfc(res0(:));
056 resWAFO = N-res+1;
057 RFMres = dtp2rfm(resWAFO,N,'CS');
058 RFMres = fliplr(RFMres)';
059 f = f-RFMres;
060 end
061
062
063
064
065 largest_res=sort([min(min(res)) max(max(res))]);
066 index_i=largest_res(1); index_j=N-largest_res(2)+1;
067 n_large_res=f(index_i,index_j); f(index_i,index_j)=0;
068 xtra_cc=[]; for i=1:n_large_res, xtra_cc=[xtra_cc largest_res]; end
069
070 r=res(:)';
071
072
073
074 undone=1;
075 for i=1:length(r)-1
076 if undone
077 if ( sum(largest_res==[r(i) r(i+1)])==2 )
078 r=[r(1:i-1) xtra_cc r(i:length(r))];
079 undone=0;
080 elseif ( sum(fliplr(largest_res)==[r(i) r(i+1)])==2 )
081 r=[r(1:i-1) fliplr(xtra_cc) r(i:length(r))];
082 undone=0;
083 end
084 end
085 end
086
087
088
089
090 if num_cc~=[]
091
092 f=f/sum(sum(f));
093
094
095
096 num_cc_limit=2*num_cc+4;
097
098
099
100 f=floor(max([1000*N 1e6])*f);
101 end
102
103
104 comb=fr2comb(f,r);
105
106
107 X=r(1);
108
109 again=1; num_runs=0;
110
111 while again
112
113 num_runs=num_runs+1;
114
115 h=0;
116
117 if r(1)>r(2)
118
119
120 if r(1)>(r(2)+1)
121 m=r(1); M=r(2); j=N-m+1; i=m-1:-1:M+1;
122 combinations=comb(i,j);
123 freq=f(i,j);
124 prob=comb2pro(combinations,freq);
125 h=min(find(rand<cumsum(prob)))-1;
126 end
127
128 if h>0
129 C=m-h;
130 i=C:m;
131 if f(C,j)>0
132 f(C,j)=f(C,j)-1;
133 i=(C+1):r(1);
134 j=N+1-r(1);
135 comb(i,j)=comb(i,j)-1;
136 r=[m-1 C r];
137 else
138 'Warning 1', break
139 r=[r(1)-1 r(2:length(r))];
140 end
141 else
142 j=N+1-r(1);
143 if r(1)>(r(2)+1)
144 i=(r(2)+1):(r(1)-1);
145 comb(i,j)=comb(i,j)-1;
146 end
147 r(1)=r(1)-1;
148
149 end
150
151 else
152
153
154 if r(1)<r(2)-1
155 m=r(2); M=r(1); i=M; j=M+1:m-1; j=N-j+1;
156 combinations=comb(i,j)';
157 freq=f(i,j)';
158 prob=comb2pro(combinations,freq);
159 h=min(find(rand<cumsum(prob)))-1;
160 end
161
162 if h>0
163 j=M:M+h; j=N-j+1;
164 c0=M+h;
165 c=N+1-c0;
166 if f(i,c)>0
167 f(i,c)=f(i,c)-1;
168 i=r(1);
169 j=N+1-((i+h-1:-1:r(1)));
170 comb(i,j)=comb(i,j)-1;
171 r=[M+1 c0 r];
172 else
173 'Warning 2', break
174 r=[r(1)+1 r(2:length(r))];
175 end
176 else
177 i=r(1);
178 j=N+1-(r(1):(r(2)-1));
179 comb(i,j)=comb(i,j)-1;
180 r(1)=r(1)+1;
181
182 end
183 end
184
185 if r(1)==r(2), r=r(2:length(r)); end
186
187 if (length(r)>=2)
188 again=1;
189 else
190 disp(' The residual is empty. Program will terminate.')
191 again=0;
192 end
193
194 X=[X r(1)];
195
196 n_tmp=length(X);
197
198
199
200 if n_tmp>2
201 if (X(n_tmp-2)<X(n_tmp-1))&(X(n_tmp-1)<X(n_tmp))
202 X=[X(1:n_tmp-2) X(n_tmp)];
203 elseif (X(n_tmp-2)>X(n_tmp-1))&(X(n_tmp-1)>X(n_tmp))
204 X=[X(1:n_tmp-2) X(n_tmp)];
205 end
206 end
207
208 comb=fliplr(triu(fliplr(comb),1));
209
210
211
212
213 if num_cc~=[]
214 if length(X)>num_cc_limit
215 disp(' The number of requested cycles has been reached.')
216 disp(' Program will terminate.')
217 again=0;
218 end
219 end
220
221 end
222
223 X = X';
224