Path: news.mathworks.com!not-for-mail
From: "Bruno Luong" <b.luong@fogale.findmycountry>
Newsgroups: comp.soft-sys.matlab
Subject: Re: interpolating/smoothing w/ monotonically increasing
Date: Thu, 13 Aug 2009 10:01:03 +0000 (UTC)
Organization: FOGALE nanotech
Lines: 183
Message-ID: <h60o8v$e66$1@fred.mathworks.com>
References: <ged003$oi7$1@fred.mathworks.com> <ged11r$f59$1@fred.mathworks.com> <ged2vl$ep1$1@fred.mathworks.com> <ged7g2$ki3$1@fred.mathworks.com> <gedajd$mp$1@fred.mathworks.com> <gedfqa$am6$1@fred.mathworks.com>
Reply-To: "Bruno Luong" <b.luong@fogale.findmycountry>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1250157663 14534 172.30.248.38 (13 Aug 2009 10:01:03 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 13 Aug 2009 10:01:03 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 390839
Xref: news.mathworks.com comp.soft-sys.matlab:562937


An anonymous reader has solicited me to post the code of monotonic interpolation. As I do not own optimization toolbox I cannot check whereas LINPROG is correctly called and provides good result. But here we go. If you see any issue with linprog call please let me know.

Bruno

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% testmonottonic.m
% Solve the problem submitted in the following thread
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/238438
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Data
% Pete Shere's (OP) test data
x=[0.1;0.2;0.3;0.4;0.5;0.6;0.7;0.8;0.9;1];

y=[183.6992;347.6811;184.4308;...
   199.1481;308.3963;453.1903;...
   490.1048;558.374;613.4362;...
   623.4727];

% L2 solution, using lsqlin, solved by John D'Errico
yl2 = [183.7
       243.75
       243.75
       243.75
        308.4
       453.19
        490.1
       558.37
       613.44
       623.47 ];
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Engine starts here, basic reshape
x = x(:); % Actually x does not play any role in the problem, it is simply
          % a correspondant abscissa to y data
y = y(:);
n = length(y);
if length(y)~=n
    error('x and y must have the same length');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% L1 fitting, BL's preference
M = eye(n);

% L1 system with slack variables u and v
%   y(i) = yfit(i) + u(i) - v(i) for all i
%   u>=0
%   v>=0
%   min sum(u+v) = sum|y-yfit|
Aeq = [M eye(n,n) -eye(n,n)];
beq = y(:);
c = [zeros(1,size(M,2)) ones(1,2*n)].';
%
LB = [-inf(1,size(M,2)) zeros(1,2*n)].';
% no upper bounds at all.
UB = inf(size(LB)); % [];

% Monotonic requirement
A = diag(ones(1,n),0)-diag(ones(1,n-1),1);
A(end,:) = [];
A(end,3*n) = 0;
b = zeros(n-1,1);

% Solve LP
% Minimizing (for x in R^n): f(x) = c'*x, subject to
%       A*x <= b        (LE)
%       Aeq*x = beq     (EQ)
%       LB <= x <= UB   (BD).
if isempty(which('linprog'))
    % Here is BL Linprog, call Matlab linprog instead to get "SOL",
    % the solution of the above LP problem
    mpsfile = 'L1_1DReg.mps';
    Contain = BuildMPS(A,b,Aeq,beq,c,LB,UB);
    OK = SaveMPS(mpsfile,Contain);
    if ~OK
        disp('Cannot write mps file');
        return
    end
    [OK outfile] = invokePCx(mpsfile);
    if ~OK
        disp('Cannot invoke PCx');
        return
    end
    [OK sol] = readPCxoutput(outfile);
    if ~OK
        disp('Cannot read PCx outfile');
        return
    end
else % Matlab linprog
    % Adjust optional LP options at your preference, see
    % http://www.mathworks.com/access/helpdesk/help/toolbox/optim/ug/linprog.html
    % http://www.mathworks.com/access/helpdesk/help/toolbox/optim/ug/optimset.html
    
    % BL has not checked because he does not have the right toolbox
    LPoption = {};
    % LPoption = {optimset(...)}    
    sol = linprog(c,A,b,Aeq,beq,LB,UB,[],LPoption{:});
end
%%% End of linprog solution

% Discard slack variables
sol = sol(1:size(M,2));
yl1 = M*sol;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Linfinity fitting
M = eye(n);

% Linfinity system with slack variables u and v
% L1 system with slack variables u and v
%   y(i)+u = yfit(i) OR
%   y(i)-v = yfit(i)
%   u>=0 and v>=0
%   min u+v = max|y-yfit|
LB = [-inf*ones(1,n) 0 0].';
UB = inf(size(LB)); % [];
A_LP=[+M -ones(n,1) zeros(n,1); ...
      -M zeros(n,1) -ones(n,1)];
b_LP=[+y;...
      -y];
f = [zeros(1,n) 1 1].';

% Monotonic requirement
A = diag(ones(1,n),0)-diag(ones(1,n-1),1);
A(end,:) = [];
A(end,n+2)=0;
b = zeros(n-1,1);

AA = [A_LP; ...
      A];
bb = [b_LP; ...
      b];
% Solve LP
% Minimizing (for x in R^n): f(x) = f'*x, subject to
%       AA*x <= bb        (LE)
%       LB <= x <= UB   (BD).
if isempty(which('linprog'))
    % Here is BL Linprog, call Matlab linprog instead to get "SOL",
    % the solution of the above LP problem    
    mpsfile='Linf_1DReg.mps';
    
    Contain=BuildMPS(AA,bb,[],[],f,LB,UB,'Linf_1DReg');
    OK=SaveMPS(mpsfile,Contain);
    if ~OK
        disp('Cannot write mps file');
        return
    end
    [OK outfile]=invokePCx(mpsfile);
    if ~OK
        disp('Cannot invoke PCx');
        return
    end
    [OK solinf]=readPCxoutput(outfile);
    if ~OK
        disp('Cannot read PCx outfile');
        return
    end
else % Matlab linprog
    % Adjust optional LP options at your preference, see 
    % http://www.mathworks.com/access/helpdesk/help/toolbox/optim/ug/linprog.html
    % http://www.mathworks.com/access/helpdesk/help/toolbox/optim/ug/optimset.html
    
    % BL has not checked because he does not have the right toolbox
    LPoption = {};
    % LPoption = {optimset(...)}
    solinf = linprog(f,AA,bb,[],[],LB,UB,[],LPoption{:});
end

yinf=solinf(1:n);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Graphic output
fig = figure(1); clf(fig);
ax = axes('Parent',fig);
plot(ax,x,y,'or',...
     x,yl1,'b',...
     x,yinf,'g', ...
     x,yl2, 'c');
legend(ax,'data','l_1','l_{inf}','l_2','Location','NorthWest');

[y yl1 yinf yl2]
%% End of the script