source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/opt/yalmip/operators/pwa_yalmip.m @ 37

Last change on this file since 37 was 37, checked in by (none), 14 years ago

Added original make3d
File size: 11.2 KB
Line 
1function varargout = pwa_yalmip(varargin)
2%PWA_YALMIP Defines a piecewise function using data from MPT
3%
4%Only intended for internal use in YALMIP
5%
6% Three classes of PWA functions can be generated
7%
8% 1) Convexity aware epigraph version with a
9%    convex domain and objective modelled as c'x <t.
10%    Note, if used in non-concvex setting, a milp
11%    model will be generated as a back-up.
12%
13%    This is the function generated to describe value
14%    function of a single mpLP problem.
15%
16% 2) Overlapping partitions with convex functions
17%    in each partition. Requires one binary per region.
18%
19%    This is the function used to describe the value
20%    function generated when not removing overlaps in
21%    binary mpLP.
22%
23% 3) Non-overlapping general pwa function. Requires
24%    one binary per region.
25%
26%    This is the function generated for value functions
27%    when remove overlaps is done. It is also used to
28%    describe the pwa optimizer.
29%
30% The input is a cell with structs. Each struct has
31% the fields Pn, Pfinal, Bi and Ci (or Fi and Gi)
32
33% Author Johan Löfberg
34% $Id: pwa_yalmip.m,v 1.3 2006/06/07 14:36:32 joloef Exp $
35
36switch class(varargin{1})
37
38    case {'struct','cell'} % Should only be called internally
39
40        if isa(varargin{2},'double')
41            % Called from YALMIP to get double
42            pwastruct = varargin{1};
43            x = varargin{2};
44            index = varargin{5};
45
46            val = inf;
47            for i = 1:length(pwastruct)
48                [ii,jj] = isinside(pwastruct{i}.Pn,x);
49                if ii
50                    for k = 1:length(jj)
51                        val = min(val,pwastruct{i}.Bi{jj(k)}(index,:)*x+pwastruct{i}.Ci{jj(k)}(index));
52                    end
53                end
54            end
55            if isinf(val)
56                val = nan;
57            end
58            varargout{1} = val;
59            return
60        end
61
62        if nargin<3
63            pwaclass = 'general'
64        end
65
66        if isa(varargin{1},'struct')
67            varargin{1} = {varargin{1}};
68        end
69
70        % Put in standard format
71        if ~isfield(varargin{1}{1},'Bi')
72            if ~isfield(varargin{1}{1},'Fi')
73                error('Wrong format on input to PWA (requires Bi or Fi)');
74            else
75                for i = 1:length(varargin{1})
76                    varargin{1}{i}.Bi = varargin{1}{i}.Fi
77                    varargin{1}{i}.Ci = varargin{1}{i}.Gi
78                end
79            end
80        end
81
82        if ~isfield(varargin{1}{1},'Pfinal')
83            error('Wrong format on input to PWA (requires field Pn)');
84        end
85
86        % Create binary variables already here to avoid generating
87        % binary variables for the same variable several times
88        switch varargin{3}
89            case 'convex'
90                % No binary needed
91                varargin{end+1} = [];
92            case 'convexoverlapping'
93                % One binary per overlap
94                varargin{end+1} = binvar(length(varargin{1}),1);
95            case 'general'
96                % one binary per region
97                varargin{end+1} = binvar(length(varargin{1}{1}.Pn),1);
98            otherwise
99        end
100
101        % Create one variable for each row
102        % Inefficient but the only way currently in YALMIP
103        varargout{1} = [];
104        varargin{end+1} = 1;
105
106        for i = 1:length(varargin{1}{1}.Ci{1})
107            varargin{end} = i;
108            varargout{1} = [varargout{1};yalmip('addextendedvariable',mfilename,varargin{:})];
109        end
110
111    case 'char' % YALMIP sends 'model' when it wants the epigraph or hypograph
112
113        switch varargin{1}
114
115            case 'graph'
116
117                % Can only generate the first class of PWA functions
118                t     = varargin{2};     % The YALMIP variables modelling this pwa
119                pwa_struct = varargin{3};% MPT structure
120                x     = varargin{4};     % Argument
121                flag  = varargin{5};     % Type of PWA function
122                d     = varargin{6};     % Binary for nonconvex cases
123                index = varargin{7};     % Which row in Bix+Ci
124
125                switch flag
126
127                    case 'convex'
128
129                        [H,K] = double(pwa_struct{1}.Pfinal);
130                        costs = reshape([pwa_struct{1}.Bi{:}]',length(x),[])'*x+reshape([pwa_struct{1}.Ci{:}]',[],1);
131                        F = set(H* x <= K) + set(costs< t);
132
133                    case 'convexoverlapping'
134
135                        % Local costs
136                        s = sdpvar(length(pwa_struct),1);
137
138                        % Some region, one defined as minimizer
139                        F = set(sum(d)==1);
140                        maxcost = -inf;
141                        mincost = inf;
142                        for i = 1:length(pwa_struct)
143                            % Reduce complexity of max function by
144                            % removing equal costs
145                            B = reshape([pwa_struct{i}.Bi{:}]',length(x),[])';
146                            C = reshape([pwa_struct{i}.Ci{:}]',[],1);
147                            [ii,jj,kk] = unique(round(1e8*[B C])/1e8,'rows');
148                            cost = reshape([pwa_struct{i}.Bi{jj}]',length(x),[])'*x+reshape([pwa_struct{i}.Ci{jj}]',[],1);
149                            [M,m] = derivebounds(cost);
150                            maxcost = max(maxcost,max(M));
151                            mincost = min(mincost,min(m));
152                            [H,K] = double(pwa_struct{i}.Pfinal);
153                            [M,m] = derivebounds(H*x - K);
154                            F = F + set(H*x-K <= (1+M)*(1-d(i)));
155                        end
156                        for i = 1:length(pwa_struct)
157
158                            B = reshape([pwa_struct{i}.Bi{:}]',length(x),[])';
159                            C = reshape([pwa_struct{i}.Ci{:}]',[],1);
160                            [ii,jj,kk] = unique(round(1e8*[B C])/1e8,'rows');
161
162                            cost = reshape([pwa_struct{i}.Bi{jj}]',length(x),[])'*x+reshape([pwa_struct{i}.Ci{jj}]',[],1);
163
164                            [M,m] = derivebounds(cost);
165                            F = F + set(cost < t + 2*(1+maxcost)*(1-d(i)));
166                        end
167                        [t_bounds] = yalmip('getbounds',getvariables(t));
168                        bounds(t,max([t_bounds(1) mincost]),min([t_bounds(2) 3*maxcost]));
169
170
171                    case 'general'
172
173                        % In one region
174                        F = set(sum(d) == 1);
175                       
176                        % Extract the wanted row
177                        for i = 1:length(pwa_struct{1}.Pn)   
178                            pwa_struct{1}.Bi{i} = pwa_struct{1}.Bi{i}(index,:);
179                            pwa_struct{1}.Ci{i} = pwa_struct{1}.Ci{i}(index);
180                        end
181                       
182                        % value in different regions
183                        costs = reshape([pwa_struct{1}.Bi{:}]',length(x),[])'*x+reshape([pwa_struct{1}.Ci{:}]',[],1);
184                        [M,m] = derivebounds(costs);
185
186                        % We know that t is always between min(m) and max(M).
187                        % However, user might know more
188                        [t_bounds] = yalmip('getbounds',getvariables(t));
189                        bounds(t,max([t_bounds(1)  min(m)]),min([t_bounds(2) max(M)]));
190
191                        % t equals some of the costs
192                        % Big-M : |costs-t| < max(costs)-min(t) < M-m
193                        F = F + set( -2*(M-m+1).*(1-d)  <= costs-t <= 2*(1+M-m).*(1-d));
194                        for i = 1:length(pwa_struct{1}.Pn)
195                            [H,K] = double(pwa_struct{1}.Pn(i));
196                            [M,m] = derivebounds(H*x - K);
197                            F = F + set(H*x - K <= M*(1-d(i)));
198                        end
199
200                    otherwise
201
202                        varargout{1} = [];
203                        varargout{2} = struct('convexity','milp','monotoncity','none','definiteness','none');
204                        varargout{3} = [];
205                        return
206                end
207
208                varargout{1} = F;
209                varargout{2} = struct('convexity','convex','monotoncity','none','definiteness','none');
210                varargout{3} = x;
211
212            case 'milp'
213
214                % FIX : Should create case for overlapping convex PWAs,
215                % used in a nonconvex fashion...
216
217                % Can only generate the first class of PWA functions
218                t     = varargin{2};     % The YALMIP variables modelling this pwa
219                pwa_struct = varargin{3};% MPT structure
220                x     = varargin{4};     % Argument
221                flag  = varargin{5};     % Type of PWA function
222                d     = varargin{6};     % Binary for nonconvex cases
223                index = varargin{7};     % Which row in Bix+Ci
224
225                switch flag
226
227                    case {'general','convex'}
228
229                        F = set(sum(d) == 1);
230                       
231                        % Extract the wanted row
232                        for i = 1:length(pwa_struct{1}.Pn)   
233                            pwa_struct{1}.Bi{i} = pwa_struct{1}.Bi{i}(index,:);
234                            pwa_struct{1}.Ci{i} = pwa_struct{1}.Ci{i}(index);
235                        end
236                       
237                        % Cost in different regions
238                        costs = reshape([pwa_struct{1}.Bi{:}]',length(x),[])'*x+reshape([pwa_struct{1}.Ci{:}]',[],1);
239                        [M,m] = derivebounds(costs);
240
241                        % We know that t is always between min(m) and max(M).
242                        % However, user might know more
243                        [t_bounds] = yalmip('getbounds',getvariables(t));
244                        bounds(t,max([t_bounds(1)  min(m)]),min([t_bounds(2) max(M)]));
245
246                        % Might be called with a convex function,
247                        % but wants the complete MILP description
248                        % Example : Someone tries to maximize value
249                        % function
250                        if length(d)~=length(costs)
251                            d = binvar(length(costs),1);
252                            F = set(sum(d) == 1);
253                        end
254                        % t equals some of the costs
255                        % Big-M : |costs-t| < max(costs)-min(t) < M-m
256                        F = F + set( -2*(M-m+1).*(1-d)  <= costs-t <= 2*(1+M-m).*(1-d));
257                        for i = 1:length(pwa_struct{1}.Pn)
258                            [H,K] = double(pwa_struct{1}.Pn(i));
259                            [M,m] = derivebounds(H*x - K);
260                            F = F + set(H*x - K <= (M+1)*(1-d(i)));
261                        end
262
263                    otherwise
264
265                        varargout{1} = [];
266                        varargout{2} = struct('convexity','convex','monotoncity','none','definiteness','none');
267                        varargout{3} = [];
268                        return
269                end
270
271                varargout{1} = F;
272                varargout{2} = struct('convexity','milp','monotoncity','milp','definiteness','milp');
273                varargout{3} = x;
274
275
276            otherwise
277                error('PWA_YALMIP called with CHAR argument?');
278        end
279    otherwise
280        error('Strange type on first argument in SDPVAR/NORM');
281end
Note: See TracBrowser for help on using the repository browser.