001 function f_rfc = mctp2rfc(f_mM,f_Mm,paramm,paramM)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022 if nargin==1
023 f_Mm=f_mM;
024 end
025
026 if nargin<3
027 paramm=[-1, 1 ,length(f_mM)];
028 paramM=paramm;
029 end
030
031 if nargin<4
032 paramM=paramm;
033 end
034
035
036 N=length(f_mM);
037 Max=sum(f_mM');
038 Min=sum(f_mM);
039 f_rfc=zeros(N,N);
040 f_rfc(N-1,1)=Max(N-1);
041 f_rfc(1,N-1)=Min(N-1);
042 for k=3:N-1
043 for i=2:k-1,
044 AA=f_mM(N-k+1:N-k+i-1,k-i+1:k-1);
045 AA1=f_Mm(N-k+1:N-k+i-1,k-i+1:k-1);
046 RAA=f_rfc(N-k+1:N-k+i-1,k-i+1:k-1);
047 nA=length(AA);
048 MA=Max(N-k+1:N-k+i-1);
049 mA=Min(k-i+1:k-1);
050 SA=sum(sum(AA));
051 SRA=sum(sum(RAA));
052 SMA=sum(MA);
053 SmA=sum(mA);
054 DRFC=SA-SRA;
055 NT=min(mA(1)-sum(RAA(:,1)),MA(1)-sum(RAA(1,:)));
056 NT=max(NT,0);
057
058 if NT>1e-6*max(MA(1),mA(1))
059
060 NN=MA-sum(AA');
061 e=(mA-sum(AA))';
062 e=flipud(e);
063 PmM=rot90(AA);
064 for j=1:nA,
065 norm=mA(nA-j+1);
066 if norm~=0
067 PmM(j,:)=PmM(j,:)/norm;
068 e(j)=e(j)/norm;
069 end
070 end
071 fx=0.;
072
073 if max(abs(e))>1e-6 & max(abs(NN))>1e-6*max(MA(1),mA(1))
074
075 PMm=AA1;
076 for j=1:nA,
077 norm=MA(j);
078 if norm~=0
079 PMm(j,:)=PMm(j,:)/norm;
080 end
081 end
082 PMm=fliplr(PMm);
083
084 A=PMm; B=PmM;
085 I=eye(size(A));
086
087 if nA==1
088 fx=NN*(A/(1-B*A)*e);
089 else
090 fx=NN*(A*((I-B*A)\e));
091 end
092 end
093
094 f_rfc(N-k+1,k-i+1)=fx+DRFC;
095
096
097
098 else
099 f_rfc(N-k+1,k-i+1)=0.;
100 end
101 end
102 m0=max(0,Min(1)-sum(f_rfc(N-k+2:N,1)));
103 M0=max(0,Max(N-k+1)-sum(f_rfc(N-k+1,2:k)));
104 f_rfc(N-k+1,1)=min(m0,M0);
105
106 end
107
108 for k=2:N
109 M0=max(0,Max(1)-sum(f_rfc(1,N-k+2:N)));
110 m0=max(0,Min(N-k+1)-sum(f_rfc(2:k,N-k+1)));
111 f_rfc(1,N-k+1)=min(m0,M0);
112 end
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139